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Abstract 

A  fundamental  requirement  in  realistic  computational  geophysical  fluid  dynamics 
is  the  optimal  estimation  of  gridded  helds  and  of  spatial-temporal  scales  directly 
from  the  spatially  irregular  and  multivariate  data  sets  that  are  collected  by  varied 
instruments  and  sampling  schemes.  In  this  work,  we  derive  and  utilize  new  schemes  for 
the  mapping  and  dynamical  inference  of  ocean  helds  in  complex  mnltiply-connected 
domains,  stndy  the  compntational  properties  of  onr  new  mapping  schemes,  and  derive 
and  investigate  new  schemes  for  adaptive  estimation  of  spatial  and  temporal  scales. 

Objective  Analysis  (OA)  is  the  statistical  estimation  of  helds  using  the  Bayesian- 
based  Gauss-Markov  theorem,  i.e.  the  npdate  step  of  the  Kalman  Filter.  The  existing 
mnlti-scale  OA  approach  of  the  Multidisciplinary  Simulation,  Estimation  and  Assim¬ 
ilation  System  consists  of  the  snccessive  ntilization  of  Kalman  npdate  steps,  one  for 
each  scale  and  for  each  correlation  across  scales.  In  the  present  work,  the  approach 
is  extended  to  held  mapping  in  complex,  multiply-connected,  coastal  regions  and 
archipelagos.  A  reasonably  accurate  correlation  function  often  reqnires  an  estimate 
of  the  distance  between  data  and  model  points,  withont  going  across  complex  land- 
forms.  New  methods  for  OA  based  on  estimating  the  length  of  optimal  shortest  sea 
paths  using  the  Level  Set  Method  (LSM)  and  Fast  Marching  Method  (EMM)  are 
derived,  implemented  and  ntilized  in  general  idealized  and  realistic  ocean  cases.  Onr 
new  methodologies  could  improve  widely-used  gridded  databases  such  as  the  climato¬ 
logical  gridded  helds  of  the  World  Ocean  Atlas  (WOA)  since  these  oceanic  maps  were 
computed  without  accounting  for  coastline  constraints.  A  new  FMM-based  method¬ 
ology  for  the  estimation  of  absolnte  velocity  nnder  geostrophic  balance  in  complicated 
domains  is  also  ontlined.  Onr  new  schemes  are  compared  with  other  approaches,  in¬ 
cluding  the  use  of  stochastically  forced  diherential  equations  (SDE).  We  hnd  that  onr 
FMM-based  scheme  for  complex,  mnltiply-connected,  coastal  regions  is  more  efficient 
and  accnrate  than  the  SDE  approach.  We  also  show  that  the  held  maps  obtained  us¬ 
ing  onr  FMM-based  scheme  do  not  reqnire  postprocessing  (smoothing)  of  helds.  The 
compntational  properties  of  the  new  mapping  schemes  are  studied  in  detail.  We  hnd 
that  higher-order  schemes  improve  the  accuracy  of  distance  estimates.  We  also  show 
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that  the  covariance  matrices  we  estimate  are  not  necessarily  positive  definite  because 
the  Weiner  Khinchin  and  Bochner  relationships  for  positive  definiteness  are  only  valid 
for  convex  simply-connected  domains.  Several  approaches  to  overcome  this  issue  are 
discussed  and  qualitatively  evaluated.  The  solutions  we  propose  include  introducing 
a  small  process  noise  or  reducing  the  covariance  matrix  based  on  the  dominant  sin¬ 
gular  value  decomposition.  We  have  also  developed  and  utilized  novel  methodologies 
for  the  adaptive  estimation  of  spatial-temporal  scales  from  irregularly  spaced  ocean 
data.  The  three  novel  methodologies  are  based  on  the  use  of  structure  functions,  short 
term  Fourier  transform  and  second  generation  wavelets.  To  our  knowledge,  this  is  the 
first  time  that  adaptive  methodologies  for  the  spatial-temporal  scale  estimation  are 
proposed.  The  ultimate  goal  of  all  these  methods  would  be  to  create  maps  of  spatial 
and  temporal  scales  that  evolve  as  new  ocean  data  are  fed  to  the  scheme.  This  would 
potentially  be  a  significant  advance  to  the  ocean  community  for  better  understanding 
and  sampling  of  ocean  processes. 

Thesis  Supervisor:  Pierre  F.  J.  Lermusiaux 
Title:  Associate  Professor 
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Chapter  1 


Introduction  and  Motivation 


The  statistical  estimation  theory  of  Objective  Analysis  (OA)  was  introduced  by 
Gandin  (1965)  to  the  held  of  meteorology  and  was  extended  to  oceanography  by 
Bretherton  et  ah  (1976).  The  theory  is  based  on  the  Gauss-Markov  theorem  (Blackett, 
1950),  and  it  provides  a  sound  basis  for  interpolating  irregularly  spaced  data  onto  a 
computational  grid.  Upto  details  of  the  set-up,  which  are  specihc  to  the  oceanic  and 
atmospheric  helds,  the  OA  scheme  is  equivalent  to  utilize  the  Kalman  update  steps 
of  the  Kalman  Filter  to  grid  the  irregularly-spaced  data.  Specihcally,  the  data  is 
gridded  based  on  the  specihed  prior  held  estimate  and  error  covariance  matrices.  The 
OA  methodology  has  been  well  formulated  for  open  oceans  without  any  landforms 
(convex  simply-connected  domains),  but  the  OA  in  complex  coastal  regions  (multiply- 
connected  domains)  is  one  of  the  ‘last’  mapping  problems  which  remains  to  be  studied 
in  detail.  This  is  one  of  the  main  research  question  of  the  present  work. 

Our  OA  research  is  carried  out  within  the  Multidisciplinary  Simulation,  Esti¬ 
mation  and  Assimilation  System  (MSEAS)  group.  MSEAS  (http://mseas.mit.edu) 
consists  of  a  set  of  mathematical  models  and  computational  methods  for  ocean  predic¬ 
tions  and  dynamical  diagnostics,  for  optimization  and  control  of  autonomous  ocean 
observation  systems,  and  for  data  assimilation  and  data-model  comparisons.  It  is 
used  for  basic  and  fundamental  research  and  for  realistic  simulations  and  predictions 
in  varied  regions  of  the  world’s  ocean,  recently  including  monitoring  (Lermusiaux, 
2007),  naval  exercises  including  real-time  acoustic-ocean  predictions  (Xu  et  ah,  2008) 
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and  environmental  management  (Cossarini  et  al.,  2009).  Several  different  models  are 
inclnded  in  the  MSEAS,  inclnding  a  new  free-snrface  primitive-eqnation  dynamical 
model  which  uses  two-way  nesting  free-snrface  and  open  boundary  condition  schemes 
(Haley  et  ah,  2008).  This  new  free-snrface  code  is  based  on  the  primitive-equation 
model  of  the  Harvard  Ocean  Prediction  System  (HOPS).  Additionally,  barotropic 
tides  are  calculated  from  an  inverse  tidal  model  (Logoutov,  2008). 

In  MSEAS,  the  Kalman  updates  for  data  gridding  are  carried  out  successively, 
from  the  largest  scale  (uniform  mean  prior)  to  the  smallest  scale,  using  a  sequential 
processing  of  observations  and  scale  separation.  In  a  two-scale  version,  a  two-staged 
OA  approach  (Lermusiaux,  1997,  1999a)  maps  the  scarcely  available  data  onto  oceanic 
helds  in  two  steps:  the  larger  and  the  smaller  scale  steps.  The  two  main  requirements 
for  the  Objective  Analysis  based  on  a  Kalman  update  (also  called  the  Gauss  Markov 
estimation  theory)  are  the  statistical  description  of  the  held  being  estimated  and  the 
observational  noise  covariance.  While  observational  noise  statistics  is  dependent  on 
the  measurement  sensor,  the  knowledge  of  the  held  statistics  does  not  come  easily 
in  oceanography  due  to  the  scarcity  of  observations.  A  description  of  held  statistics 
is  often  provided  by  a  simple  analytical  correlation  function  which  depends  on  the 
spatial  separation  distance  and  the  spatial-temporal  scales  (Carter  and  Robinson, 
1987).  Other  schemes  also  utilize  dynamical  models  to  construct  covariances. 

Our  research  study  on  Objective  Analysis  for  coastal  regions  has  been  motivated 
by  the  Philippines  Straits  Dynamics  Experiment  (PhilEx)  sponsored  by  the  Office  of 
Naval  Research.  The  goal  of  PhilEx  is  to  enhance  understanding  of  the  oceanographic 
processes  and  features  arising  in  and  around  straits,  and  to  improve  the  capability 
to  predict  the  inherent  spatial  and  temporal  variability  of  these  regions  using  models 
and  advanced  data  assimilation  techniques.  There  are  several  examples  of  Objective 
Analysis  in  coastal  regions  (Hessler,  1984,  Stacey  et  ah,  1988,  Paris  et  al.,  2002),  but 
the  methodologies  employed  in  these  examples  do  not  satisfy  coastline  constraints 
(e.g.  there  should  be  no  direct  relationship  across  landforms). 

New  methodologies  for  held  (e.g.  temperature,  salinity,  biology,  and  velocity) 
mapping  in  complex  multiply-connected  coastal  domains  and  archipelagos  are  derived 
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and  demonstrated  in  this  work.  These  methodologies  will  likely  be  very  useful  in  im¬ 
proving  the  World  Ocean  Atlas  (WOA)  climatologies  in  complex  multiply-connected 
domains.  The  WOA  provides  global  ocean  climatology  containing  monthly,  seasonal 
and  annual  means  of  temperature  (T)  and  salinity  (S)  helds  at  standard  ocean  depths. 
The  temperature  and  salinity  climatologies  presented  as  part  of  the  WOA  (Levitus, 
1982),  which  is  also  termed  as  ‘Levitus  Climatology’  and  its  atlas  updates  in  1994 
(Levitus  and  Boyer,  1994,  Levitus  et  ah,  1994),  1998  (Antonov  et  ah,  1998a,b,c,  Boyer 
et  al.,  1998a,b,c),  2001  (Stephens  et  ah,  2002,  Boyer  et  ah,  2002)  and  2005  (Locarnini 
et  ah,  2006,  Antonov  et  ah,  2006,  Garcia  et  al.,  2006a,b)  have  proven  to  be  valuable 
tools  for  studying  the  hydrographic  structures  of  the  World’s  oceans.  The  WOA 
climatologies  have  been  very  useful  for  providing  initial  and  boundary  conditions  to 
ocean  circulation  models.  As  its  MSEAS  counterpart,  the  OA  procedure  for  ‘Levi¬ 
tus  Climatology’  requires  the  use  of  an  analytical  correlation  function  to  determine 
the  covariance  (or  weight  function,  as  described  by  Levitus  (1982)).  If  the  “straight 
Euclidean  distance”  (the  straight  line  distance  between  two  points)  is  used  in  such 
analytical  correlation  functions,  the  distance  estimate  is  inappropriate  for  complex 
multiply-connected  domains,  as  this  “straight  Euclidean  distance”  goes  across  land 
and  so  violate  the  coastline  constraints.  An  appropriate  measure  of  distance  to  be 
used  in  the  correlation  function  for  OA  in  such  complex  multiply-connected  regions 
should  be  longer.  It  is  nonetheless  the  length  of  the  optimal  shortest  sea  path  i.e., 
the  shortest  path  without  going  across  complex  landforms. 

Such  an  optimal  shortest  sea  path  in  complex  multiply-connected  regions  can 
be  obtained  using  the  following  numerical  techniques:  the  Level  set  method  (LSM) 
(Osher  and  Sethian,  1988,  Sethian,  1999b)  and  the  Fast  Marching  Method  (EMM) 
(Sethian,  1996,  1999b).  These  methods  model  the  propagation  of  evolving  boundaries 
using  appropriate  PDF’s.  They  have  been  applied  in  both  the  Philippines  Archipelago 
and  Dabob  Bay  (WA,  USA)  regions.  They  are  also  compared  to  the  SDE  approach 
proposed  by  Lynch  and  McGillicuddy  (2001).  Other  optimization  methods  for  path 
planning,  for  example  Dikjstra’s  algorithm  (Bertsimas  and  Tsitsiklis,  1997)  and  Bre- 
senham  line  algorithm  (Bresenham,  1965)  could  also  be  used  for  mapping  in  complex 
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domains,  but  our  FMM  and  LSM  schemes  are  shown  to  be  computationally  more 
efficient  and  more  accurate.  The  optimal  path  length  can  also  be  used  in  conjunction 
with  the  methodology  proposed  by  Lermusiaux  et  ah  for  three-dimensional,  mul¬ 
tivariate  and  multi-scale  spatial  mapping  of  geophysical  helds  and  their  dominant 
errors  (Lermusiaux  et  ah,  1998,  2000,  Lermusiaux,  2002)  in  coastal  regions.  This 
method  reduces  the  dimension  of  the  error  covariance  matrices  by  focusing  on  the 
error  subspace  formed  by  dominant  eigen-decomposition  of  the  a-priori  covariance 
(Lermusiaux  and  Robinson,  1999). 

The  new  LSM  and  FMM  schemes  are  also  used  in  this  work  to  estimate  the  mini¬ 
mum  vertical  area  along  any  path  between  two  islands.  Vertical  areas  across  landforms 
are  needed  to  compute  the  transport  streamfunction  along  the  island  coastlines,  which 
minimizes  the  inter  island  transport.  Such  estimates  of  the  transport  streamfunction 
will  aid  in  the  computation  of  absolute  velocity  under  geostrophic  balance  (Wunsch, 
1996)  in  complex  domains  with  islands. 

Computational  properties  of  the  new  mapping  schemes  are  also  investigated  in  de¬ 
tail.  To  reduce  the  computational  cost  and  to  understand  the  impact  of  the  individual 
data,  sequential  processing  of  observations  (Parrish  and  Cohn,  1985,  Cho  et  ah,  1996) 
is  discussed  utilized.  By  dehnition,  the  prior  covariance  matrix  should  be  positive 
dehnite.  According  to  the  Wiener-Khinchin  and  Bochner  theorem  (Papoulis,  1991, 
Yaglom,  2004,  Dolloff  et  ah,  2006),  the  covariance  matrix  based  on  analytical  correla¬ 
tion  function  will  be  positive  dehnite  if  a  Fourier  transform  (or  the  spectral  density  of 
the  correlation  function)  is  non-negative  for  all  frequencies.  These  theorems  are  valid 
only  for  convex  simply-connected  domains.  In  our  complex  multiply-connected  do¬ 
mains,  the  covariance  matrix  may  become  negative  due  to:  a.  Numerical  error  in  the 
computation  of  the  optimal  path  length  using  our  new  FMM/LSM  based  schemes  b. 
The  presence  of  landforms.  These  issues  may  lead  to  divergence  problems  (Brown  and 
Hwang,  1997)  of  held  mapping  schemes.  Therefore,  the  following  two  questions  were 
resolved  and  investigated:  a).  What  are  the  computational  errors  in  optimal  path 
lengths  computed  using  the  FMM/LSM  and  how  can  they  be  reduced?  b).  What  are 
the  computational  issues  including  non-positive  dehnite  covariances  that  arise  in  a 
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multiply-connected  coastal  domain  and  how  can  they  be  remedied?  These  computa¬ 
tional  studies  were  indispensable  for  the  development  of  our  novel  FMM/LSM  based 
scheme  for  complex  multiply-connected  domains. 

Apart  from  the  OA  research,  another  major  question  for  ocean  studies  is  the 
estimation  of  spatial  and  temporal  scales  from  irregular  ocean  held  measurements. 
Such  scale  estimation  is  potentially  a  signihcant  advance  to  the  ocean  community 
for  better  understanding  and  sampling  of  ocean  processes.  The  issue  is  that  this  is  a 
challenging  problem,  especially  because  of  the  irregular  properties  of  the  data  but  also 
due  to  the  multi-scale  turbulent  and/or  intermittent  ocean  dynamics.  The  knowledge 
of  the  spatial-temporal  scales  is  linked  to  our  OA  results:  it  provides  a  direct  estimate 
of  parameters  needed  for  the  analytical  correlation  function.  These  scales  can  also 
be  used  to  obtain  a  measure  of  the  internal  Rossby  radius  of  deformation,  another 
relevant  length  scale  in  atmospheric  and  ocean  sciences.  The  internal  Rossby  radius  is 
the  length  scale  at  which  rotational  effects  become  as  important  as  buoyancy  effects. 

We  have  proposed  and  implemented  novel  methods  for  adaptive  spatial-temporal 
scale  estimation  from  irregular  held  measurements.  The  ultimate  goal  of  these  new 
methods  would  be  to  create  maps  of  spatial  and  temporal  scales  that  evolve  as  ocean 
data  are  collected  or  are  fed  to  to  the  scale  estimation  scheme,  all  without  having 
to  map  the  data.  Denman  and  Freeland  (1985)  proposed  an  approach  for  estimating 
scales  using  the  so-called  structure  function.  We  have  further  developed  this  approach, 
which  is  based  on  using  the  structure  function  and  the  non-linear  least  square  htting 
methods,  to  obtain  adaptive  scale  estimates  that  vary  in  space  and  time.  If  the  signal 
is  stationary,  i.e.  the  scales  do  not  change  with  time,  then  the  Fourier  analysis  can 
be  very  useful  tool  for  scale  estimation.  Fourier  analysis  is  valid  only  for  stationary 
signals,  but  short  term  Fourier  transforms  (STFT)  can  be  used  for  non-stationary 
signals  as  well,  even  though  it  has  some  resolution  problems.  We  have  proposed  a 
new  method  for  adaptive  scale  estimation  based  on  STFT  and  have  illustrated  the 
method  using  chirp  signal  data.  Wavelet  Analysis  is  another  approach  which  helps  in 
overcoming  the  resolution  issues,  and  which  can  potentially  be  very  useful  for  adaptive 
scale  estimation.  However,  the  limitation  of  classic  wavelet  analysis  and  STFT  is  in 
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dealing  with  the  irregularly  sampled  ocean  data.  Using  adaptive  methods  based  on 
second  generation  wavelets  (Sweldens,  1998,  Jansen  and  Oonincx,  2005),  which  are 
applicable  to  irregularly  sampled  and  non-dyadic  data  sets,  to  learn  the  largest  and 
the  most  energetic  scales  from  the  irregularly  spaced  ocean  data  is  discussed  and 
exemplihed  using  a  chirp  signal  data  set. 

In  Chapter  2,  we  review  the  two  staged  multi-scale  static  held  mapping  approach 
from  MSEAS  and  the  Objective  Analysis  scheme  used  for  ‘Levitus  Climatology’. 
In  Chapter  3,  we  derive  and  utilize  the  new  OA  methodologies  based  on  the  Level 
Set  Method  and  the  Fast  Marching  Method  and  compare  them  to  existing  schemes 
including  the  OA  based  on  the  stochastically  forced  differential  equations.  A  new 
optimization  approach  is  proposed  for  computing  the  transport  streamfunction  and 
absolute  velocity  under  geostrophic  balance  by  minimizing  the  inter-island  transport. 
In  Chapter  4,  applications  of  our  new  methodologies,  for  the  complex  regions  of 
Dabob  Bay  and  Philippines  Archipelago  are  presented.  In  Chapter  5,  we  study  the 
computational  properties  of  our  new  mapping  schemes.  Chapter  6  introduces  our 
novel  methodologies  for  adaptive  scale  estimation  including  the  use  of  the  structure 
function,  the  use  of  STFT  and  the  use  of  second  generation  wavelets.  Chapter  7 
consists  of  a  summary  and  conclusions. 
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Chapter  2 


Objective  Analysis  Approach 


The  multi-scale  OA  approach  from  MSEAS  (Section  2.1)  and  the  approach  for 
‘Levitus  Climatology’  (Section  2.2)  for  held  mapping  are  summarized  in  this  Chapter. 
These  approaches,  which  require  the  computation  of  Euclidean  distance,  are  well 
established  for  mapping  heterogeneous,  multivariate,  irregular  data  (Gandin,  1965, 
Bretherton  et  ah,  1976,  Carter  and  Robinson,  1987,  Daley,  1993)  in  open  oceans 
without  islands  or  archipelagos  as  well  as  in  atmospheric  sciences. 


2.1  Multi-scale  static  field  estimation:  MSEAS 
Objective  Analysis 

The  two  staged  OA  approach  (Lermusiaux,  1997,  1999a)  in  MSEAS,  utilizes  the 
Gauss-Mar kov  or  minimum  error  variance  criterion  (Blackett,  1950)  to  map  observa¬ 
tions  to  the  numerical  grid.  Let  us  denote  the  vector  of  numerical  grid  point  locations 
as  X  and  the  vector  of  measurement  locations  as  X,  then  the  OA  estimate  of  the  held 
(.^OA)  Q]2  background  held  ('^,d)  is  given  by: 

_  ^  +  C'or(x,X)[C'or(X,X)  +  R]-i[d-d] 

=  ^  +  K[d-d]  (2.1) 
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where,  d  =  H  is  the  observation  matrix,  d  is  the  sensor  data  vector,  R  is  the 
observational  error  covariance  matrix,  and  Gain  (K)  is  given  by: 

K  =  Gor(x,X)[C'or(X,X) +  R]-^  (2.2) 

The  error  covariance  of  the  estimated  held  is  given  by: 

=  E[(x-i?[x])(x-i?[x])^] 

=  Gor(x,  x)  —  Gor(x,  X)[Gor(X,  X)  +  R]“^Gor(X,  x) 

=  Gor(x,  x)  —  KGor(X,  x).  (2.3) 

A  comparison  between  the  MSEAS  update  equations  (OA)  and  the  Kalman  hlter 
update  equations  is  made  in  Table  2.1. 


Kalman  Filter  Update  Equations 

MSEAS  Update  equations  (OA) 

Kalman  gain: 

OA  estimator: 

Gain  =  C'or(x,  X)[C'or(X,X)  +  R]”^ 

State  estimate  update  equation: 

xt  =  xt\t-i  +  Kt{yt  -  HtXtit-i) 

State  estimate  update  equation: 

=  '{p  Gain[d  —  d]. 

Error  Covariance  equation: 

Pt  =  il- 

Error  Covariance  equation: 

pOA  ^  Cor(x,  x)  —  Gain  x  Cor(X,x) 

Table  2.1:  Comparison  between  the  Kalman  Filter  and  MSEAS  OA  update  equations 


Thus,  the  update  equations  of  OA  are  equivalent  to  the  update  equations  of  the 
discrete  Kalman  hlter  (KF)  algorithm  where  the  background  error  correlation  matrix 
for  the  held-to-data  points,  Cor(x,X),  and  the  background  correlation  matrix  at  the 
data  points,  Cor(X,X),  are  directly  related  to  the  KF  a  priori  error  covariance  matrix 
Pt\t-i  i-e.  Cor(x,X)  =  and  Cor(X,X)  =  HtPt\t-iHj  {Ht  is  the  observation 

matrix).  The  matrix  R  represents  the  error  covariance  for  the  sensor  data  d  at  the 
data  points.  This  matrix  is  often  chosen  diagonal  with  a  uniform  non-dimensional  er¬ 
ror  variance  e^,  i.e.  R  =  e^J.  In  MSEAS,  the  correlation  matrices  are  often  generated 
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from  the  following  isotropic  function: 


Cor{r) 


/  ^2 

■  ^  ,r2  AtA' 

1  exp 

0.5  X  (^2  +  ^2  ) 

(2.4) 


Here,  At  is  the  difference  between  the  observation  and  the  estimation  time  and  r  is 
the  decorrelation  time  scale.  The  parameters  Lq  and  Lg  are  the  zero-crossing  and  the 
e-folding  length  scales.  The  scalar  r  is  the  spatial  separation  distance.  The  Euclidean 
distance  does  not  satisfy  coastline  constraints.  So,  the  use  of  the  optimal  distance 
(the  minimum  distance  between  two  points  without  going  across  complex  landforms) 
is  proposed.  LSM  or  EMM  can  be  utilized  to  obtain  such  optimal  distance  between 
any  two  points  in  a  complex  (for  example,  multi-island)  multiply-connected  coastal 
region.  These  methods  are  efficient  and  accurate  solvers  for  optimal  distances  in 
multiply-connected  complex  domains  and  satisfy  coastline  constraints. 

The  MSEAS  OA  is  carried  out  in  two  stages.  In  the  hrst  stage,  the  largest  dynam¬ 
ical  scales  are  mapped  onto  computational  grid  using  the  parameters  (r,  Lq,  Te)LS- 
The  background  held  for  this  stage  is  often  chosen  to  be  constant  and  equal  to  the 
horizontal  mean  of  all  the  observations.  In  the  second  stage,  the  smaller  scales  are 
mapped  using  the  coefficients  (r,  Lq,  Le)ME  often  corresponding  to  the  most  energetic 
(meso)  scales.  The  background  held  for  this  stage  is  the  hrst  stage  OA.  The  block 
diagram  for  the  two  staged  MSEAS  OA  is  shown  in  Figure  C-1.  A  major  assumption 
in  this  OA  approach  is  that  the  errors  in  the  largest  and  the  most  energetic  stages  are 
statistically  independent.  The  accuracy  of  the  held  estimates  obtained  using  OA  also 
depends  on  the  spatial  and  time  scale  parameters  used  in  the  analytical  correlation 
function,  as  well  as  the  correlation  function  itself. 


2.2  Objective  Analysis  approach  for  the  ‘Levitus 
Climatology’ 

The  objective  analysis  scheme  used  for  ‘Levitus  Climatology’  (Levitus,  1982,  Lo- 
carnini  et  ah,  2006,  Antonov  et  ah,  2006,  Garcia  et  ah,  2006a,b)  has  its  origins  in  the 
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work  of  Cressman  (1959)  and  Barnes  (1964).  This  scheme  is  based  on  adding  “correc¬ 
tions”,  which  are  computed  as  a  distance- weighted  mean  of  all  grid  point  difference 
values,  to  the  first-guess  held.  Initially,  the  World  Ocean  Atlas  1994  (WOA94)  used 
the  Barnes  (1973)  scheme  which  requires  only  a  single  “correction”  to  the  hrst-guess 
held  at  each  grid  point  in  comparison  to  the  successive  correction  method  of  Cress- 
man  (1959)  and  Barnes  (1964).  This  was  done  to  reduce  the  computing  time.  Barnes 
(1994)  suggests  using  the  multi-pass  analysis  when  computing  time  is  not  an  issue. 
The  analysis  scheme  used  in  WOA98,  WOAOl  and  WOA05  is  a  three-pass  “correc¬ 
tion”  scheme.  The  inputs  to  this  analysis  scheme  are  one-degree  square  means  of  the 
observed  data  values,  and  a  hrst-guess  held.  The  diherence  between  the  observed 
mean  and  the  hrst-guess  held  is  then  computed.  An  inhuence  radius  is  specihed 
next  and  a  correction  to  the  hrst-guess  value  at  all  the  grid  points  is  computed  as  a 
distance-weighted  mean  of  all  the  grid  point  diherence  values  that  lie  within  the  area 
around  the  grid  point  dehned  by  the  inhuence  radius.  Mathematically,  the  correction 
factor  derived  by  Barnes  (1964)  is  given  by: 


a 


Eti  H-.Q, 
Ef.,  1-14 


(2.5) 


where, 

(i,j)  -  coordinates  of  a  grid  point  in  east-west  and  north-south  directions  respectively; 
Cij  -  correction  factor  at  the  grid  point  coordinates 

d  -  the  number  of  data  points  that  fall  within  the  area  around  point  {i,j)  dehned  by 
the  inhuence  radius; 

Qs  -  diherence  between  the  observed  mean  and  the  hrst-guess  at  the  data  point 
in  the  inhuence  area; 

Ws  =  exp{—Er‘^/R‘^)  (for  r  <  R]  Wg  =  0  for  r  >  R)  =  Correlation  weight; 
r  -  distance  of  observation  from  grid  point; 

R  -  inhuence  radius; 

E  =  A. 

At  each  grid  point,  the  analyzed  value  Gij  is  the  sum  of  the  hrst  guess  Eij  and 


the  correction  Qj.  The  expression  is: 


Gij  —  Fij  +  Cij  (2-6) 

If  there  is  no  data  within  the  area  dehned  by  the  influence  radius,  the  correction  is 
zero  and  the  analyzed  value  of  the  held  is  the  same  as  the  hrst-guess.  The  analysis 
scheme  is  set  up  such  that  the  inference  radius  can  be  varied  in  each  iteration.  To 
progressively  analyze  the  smaller  scale  phenomena  with  each  iteration,  the  analysis 
begins  with  a  large  inference  radius  which  is  decreased  gradually  with  each  iteration. 

Equation  2.6  can  also  be  expressed  in  the  matrix  form,  which  is  given  by 

G  =  F  +  [dmc/(Werf)]-^WQ  (2.7) 

Here  n  is  the  number  of  model  points,  the  analyzed  held  G  and  the  hrst  guess  F  are 
n-by-1  vectors,  the  correlation  weight  matrix  W  is  a  n-by-d  matrix,  the  diherence 
between  the  observed  mean  and  the  hrst-guess  at  the  data  point  Q  is  a  d-by-1  vector 
and  Gd  is  a  d-by-1  vector  with  unit  entities.  The  operation  diag(v)  creates  a  diagonal 
matrix  i.e.  it  puts  the  vector  v  on  the  main  diagonal. 

Analogous  to  the  Kalman  Gain  (K)  from  the  Gauss  Markov  criterion  {K  = 
C'or(x,  X)[C'or(X,X)  -|-  R]“^),  Equations  2.7  and  2.1  show  that  a  similar  Gain  ma¬ 
trix  {Kl  =  [(im5f(Werf)]“^W)  can  be  dehned  for  the  Levitus  methodology.  While  the 
multi-scale  OA  approach  in  MSEAS  is  based  on  Gauss  Markov  estimation  theory, 
the  Levitus  OA  is  based  on  estimating  the  held  by  computing  the  distance-weighted 
mean  of  all  grid  point  diherence  values  (between  the  mean  and  hrst-guess  held)  in 
the  inference  radius  and  then  adding  it  to  the  hrst-guess  held.  Thus,  the  choice  of 
hrst  guess-held  is  very  important  in  the  ‘Levitus  OA’  analysis.  On  the  other  hand, 
in  Gauss  Markov  estimation,  the  hrst-guess  held  is  often  the  mean  of  the  data  values 
and  the  correction  is  made  in  the  Kalman  update  step  by  computing  the  diherence 
between  the  data  and  the  interpolated  value  of  the  hrst-guess  on  the  data  location. 
The  Gauss  Markov  estimation  theory  also  requires  the  knowledge  of  error  covariance 
of  the  observation  noise  (R). 
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The  above  approaches  for  MSEAS  OA  and  ‘Levitus  Climatology’,  which  are  based 
on  computing  the  covariance  or  the  weight  factors  by  providing  Euclidean  distance  as 
an  input  to  the  correlation  function,  are  valid  only  for  open  oceans.  New  methodolo¬ 
gies  based  on  computing  the  length  of  the  optimal  path  using  the  Level  Set  Method 
and  the  Fast  Marching  Method  are  discussed  in  Chapter  3. 
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Chapter  3 


OA  methodologies  in  complex 
coastal  regions  and  archipelagos 


This  chapter  describes  the  new  methodologies  for  Objective  Analysis  in  complex 
multiply-connected  coastal  regions  and  compares  them  to  previous  methods.  The 
new  methodologies  are  based  on  hnding  the  length  of  the  optimal  path  using  the 
Level  Set  Method  and  the  Fast  Marching  Method,  as  described  in  Section  3.1.  In 
Section  3.2,  we  discuss  other  methods  including  the  methodology  proposed  by  Lynch 
and  McGillicuddy  (2001)  based  on  using  the  stochastically  forced  differential  equation 
(SDE)  for  the  held  and  the  methodology  based  on  using  the  SDE  for  the  covariance 
(Logoutov,  Lermusiaux,  personal  communication).  An  improved  version  of  the  op¬ 
timization  methodology  proposed  by  Haley  and  Lermusiaux  (2009),  for  estimating 
the  inter-island  transport  and  for  estimating  the  absolute  velocity  under  geostrophic 
balance  in  multiply-connected  coastal  regions  is  presented  in  Section  3.3. 


3.1  Novel  Objective  Analysis  Methodologies  based 
on  optimal  path  length 

The  optimal  path  in  a  domain  having  complex  multiply-connected  islands  and 
archipelagos  is  dehned  as  the  shortest  path  between  two  points  without  going  across 
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landforms.  Optimization  method  like  Dijkstra’s  algorithm  (Bertsimas  and  Tsitsiklis, 
1997)  can  be  utilized  to  obtain  the  optimal  path  length.  Apart  from  being  computa¬ 
tionally  expensive,  Dijkstra’s  method  is  also  inaccurate.  The  accuracy  of  FMM  and 
LSM  is  illustrated  in  Figure  C-2  where  it  is  compared  to  the  optimization  method 
(Dijkstra’s  method  with  norm  p  =  2)  for  optimal  distance  computations  in  a  complex 
domain  (Takei,  2006).  Another  classic  method  for  approximately  computing  the  opti¬ 
mal  path  length  is  the  Bresenham  line  algorithm  (Bresenham,  1965).  This  algorithm 
is  used  for  computer  control  of  a  digital  plotter.  In  the  oceanic  context,  the  modified 
Bresenham  line  algorithm  can  be  used  to  obtain  the  length  of  the  optimal  path.  The 
modihed  algorithm  is  explained  using  Figure  C-3,  which  has  a  complex  island  in  the 
domain.  To  compute  the  length  of  the  optimal  path  between  the  points  A  and  B  in 
Figure  C-3,  a  straight  line  connecting  the  two  points  is  drawn.  Since  the  optimal  path 
does  not  go  across  landforms,  the  straight  line  path  is  modihed  such  that  it  goes  along 
the  boundary  of  the  island  upon  hitting  the  island.  The  optimal  path  between  the 
points  A  and  B  is  drawn  in  Figure  C-3.  The  limitation  of  the  Bresenham  algorithm, 
apart  from  being  computationally  expensive,  is  that  the  optimal  path  length  com¬ 
puted  using  this  algorithm  is  discontinuous.  Again,  this  is  illustrated  in  Figure  C-3. 
The  straight  line  between  points  A  and  C  does  not  intersect  the  complex  island,  and 
therefore  the  straight  line  path  is  the  optimal  path.  However,  it  is  clearly  observed 
from  Figure  C-3  that  the  length  of  the  optimal  path,  which  is  computed  from  the 
Bresenham  line  algorithm,  between  points  A  and  B  is  signihcantly  larger  than  the 
length  of  the  optimal  path  between  points  A  and  C.  Since,  points  B  and  C  are  not 
very  far  from  each  other,  the  length  of  the  true  optimal  path  between  points  A  and 
B  and  between  points  A  and  C  should  not  differ  signihcantly.  This  example  clearly 
illustrates  the  limitation  of  the  Bresenham  line  algorithm  for  computing  the  optimal 
path  in  the  oceanic  context. 

We  propose  to  utilize  the  methodology  based  on  using  level  sets  for  computing 
the  length  of  the  optimal  path.  The  novel  Objective  Analysis  methodology  based  on 
computing  the  length  of  the  optimal  path  using  the  Level  Set  Method  is  discussed 
in  Section  3.1.1  and  the  methodology  based  on  using  the  Fast  Marching  Method  is 
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discussed  in  Section  3.1.2.  These  methods  can  efficiently  and  accurately  compute  the 
optimal  path  length. 


3.1.1  Objective  Analysis  using  the  Level  Set  Method  (LSM) 

A  level  set  of  a  real- valued  function  0  of  n  variables  is  a  set  of  the  form: 

{(Xi,  ...,Xn)  =  C}  (3.1) 

where,  c  is  a  constant.  That  is,  a  level  set  is  the  set  where  the  function  0  takes  on  a 
given  constant  value  c. 

Osher  and  Sethian  (1988),  Sethian  (1999b)  proposed  a  numerical  technique,  which 
is  called  the  Level  Set  Method,  to  implicitly  represent  and  model  the  propagation  of 
evolving  interfaces  under  the  influence  of  a  given  velocity  held  using  appropriate 
partial  differential  equations  (PDE’s).  An  initial  value  formulation  describing  the 
interface  motion  is  now  discussed.  The  initial  position  of  interfaces  are  given  by  level 
sets  of  the  function  0.  The  evolution  of  this  function  0  is  linked  to  the  propagation 
of  the  interface  through  a  time-dependent  level  set  equation.  Interfaces  can  be  rep¬ 
resented  explicitly  (parametrized  interfaces  i.e.  interfaces  given  by  x  =  x(s),  where 
s  is  the  parameter)  or  implicitly  (interfaces  given  by  zero  level  set  i.e.  0(x)  =  0). 
Using  the  implicit  representation  0(x),  where  x  is  the  position  vector,  the  convection 
equation  can  be  solved  to  propagate  level  sets  by  a  velocity  held  v; 

0t  -I-  V.V0  =  0  (3.2) 

In  many  cases,  one  is  interested  only  in  the  motion  normal  to  the  boundary.  Therefore, 
the  velocity  v  can  be  represented  using  the  scalar  speed  function  F  and  the  normal 
direction  n.  Thus. 


V  =  Fn 


pit 

|V0| 


(3.3) 
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The  hyperbolic,  non-linear  (Hamilton-Jacobi  equation)  level  set  equation,  obtained 
from  Equations  3.2  and  3.3,  is  given  by: 

0t  +  F|V0|=O  (3.4) 

The  following  Erst  order  upwinded  finite  difference  approximation  can  be  used  to 
numerically  solve  the  level  set  equation  (2-dimensional  in  space)  (Osher  and  Sethian, 
1988,  Sethian,  1999b): 


where. 


At[max{F,  +  min{F,  0)Vjj 


V+-  =  0)^  -h  0)^ 

max{D~^ (f)F^  0)^  -|-  min{D~^^(l)F,  0)^]^'^^ 

=  [min{D~'^(l)F,  0)^  -h  max{D^'^(()F,  0)^  -h 

min{D~'^(j)F^  0)^  -|-  max{D^'^(j)'^  -,  (3.5) 


Here,  D“^  is  the  first  order  backward  difference  operator  in  the  x-direction;  is 
the  first  order  forward  difference  operator  in  x-direction,  etc.  Mathematically,  these 
operators  are  given  by: 


D- 


Ax 


D 


+x^ 


Ax 


(3.6) 


The  level  set  equation  is  an  initial  value  problem  which  tracks  the  evolution  of  the 
level  sets  (;/)=constant  assuming  F  is  given  by  the  specifics  of  the  evolution  of  the  0 
for  a  particular  problem. 


If  the  scalar  speed  function  of  the  front  F  is  non-negative,  then  the  steady  state 
boundary  value  problem,  known  as  the  Eikonal  equation,  can  be  formulated  to  evalu¬ 
ate  the  arrival  time  function  T(x).  The  Eikonal  equation  representing  the  time  T(x) 
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for  the  “frontal  interface”  to  reach  the  position  x  from  its  initial  position  is  given  by: 


F|VT|  =  1  (3.7) 

The  Eikonal  equation  simply  states  that  the  gradient  of  the  arrival  time  function  is 
inversely  proportional  to  the  speed  of  the  front.  To  solve  the  Eikonal  equation,  a 
time  dependent  problem  is  proposed.  The  time  evolved  steady  state  solution  of  the 
resultant  Hamilton-Jacobi  equation  is  the  Eikonal  equation.  Mathematically,  this  is 
written  as: 


Tt  +  F|  VT|  =  1  F|VT|  =  1  (3.8) 

This  Hamilton-Jacobi  equation  (Equation  3.8  (Left))  can  be  discretized  using  the  nu¬ 
merical  scheme  for  the  Level  Set  equation.  The  steady  state  solution  of  this  Hamilton- 
Jacobi  equation  will  be  the  solution  of  the  Eikonal  equation  (Equation  3.8  (Right)). 

The  Level  Set  Method  has  been  used  in  a  wide  variety  of  applications  which  in¬ 
clude  the  arrival  time  problems  in  the  control  theory,  generation  of  minimal  surfaces, 
flame  propagation,  fluid  interfaces,  shape  reconstruction  etc.  In  the  oceanic  context, 
the  method  can  be  used  to  determine  the  optimal  distance  dehned  here  as  the  mini¬ 
mum  distance  between  two  points  without  going  across  complex  landforms. 
Numerics  and  operation  connt  for  the  LSM:  MATLAB  code  has  been  devel¬ 
oped  for  Objective  Analysis  using  the  Level  Set  Method.  For  estimating  the  optimal 
distance,  the  scalar  speed  function  F  is  set  to  0  for  the  grid  points  on  the  land  and  1 
for  the  grid  points  on  the  water.  The  level  set  T(x),  which  is  the  arrival  time  func¬ 
tion,  also  represents  the  optimal  distance  from  the  starting  position  to  the  position 
vector  X  for  the  above  speed  function  F .  The  above  OA  approach,  which  is  based  on 
computing  the  evolution  of  all  the  level  sets  and  not  simply  the  zero  level  set  corre¬ 
sponding  to  the  front  itself,  has  an  operation  count  of  O(N^)  in  two  dimensions  for 
grid  points  (Sethian,  1999b).  Thus,  it  is  a  computationally  expensive  technique 
since  an  extra  dimension  has  been  added  to  the  problem. 

A  modihed  approach  named  ‘Fast  Marching  level  set  method’,  which  signihcantly 
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reduces  the  operation  count,  is  described  in  the  next  section.  Roughly  speaking,  the 
two  possible  ways  to  view  these  solution  techniques  are  either  iteration  towards  the 
solution,  or  direct  construction  of  the  stationary  solution  T.  While  LSM  constructs 
the  solution  to  the  Eikonal  equation  (Equation  3.7)  by  iterating  towards  the  solution, 
EMM  is  based  on  direct  construction  of  the  stationary  solution  T. 

3.1.2  Objective  Analysis  using  the  Fast  Marching  Method 
(FMM) 

The  Fast  Marching  Method  (FMM)  for  monotonically  advancing  fronts,  which  has 
been  proposed  by  Sethian  (1996,  1999b),  is  described  in  this  section.  This  method 
leads  to  an  extremely  fast  scheme  for  solving  the  Eikonal  equation  (Equation  3.7). 
The  Level  set  method,  which  is  described  in  Section  3.1.1,  relies  on  computing  the 
evolution  of  all  level  sets  by  solving  an  initial  value  partial  differential  equation  us¬ 
ing  numerical  techniques  from  hyperbolic  conservation  laws.  As  an  alternative,  an 
efficient  modihcation  is  to  perform  the  work  only  in  the  neighborhood  of  the  zero 
level  set,  as  this  is  known  as  the  ‘narrow  band  approach’.  The  basic  idea  of  this 
alternative  approach  is  to  tag  the  grid  points  as  either  “alive” ,  “land  mines”  or  “far 
away”  depending  on  whether  they  are  inside  the  band,  near  its  boundary,  or  outside 
the  band,  respectively.  The  work  is  performed  only  on  alive  points,  and  the  band  is 
reconstructed  once  the  land  mine  points  are  reached. 

FMM,  which  allows  boundary  value  problems  to  be  solved  without  iterations,  is 
now  discussed  in  detail.  The  method  is  applicable  to  monotonically  advancing  fronts 
(i.e.  the  front  speed  (F  >  0  or  F  <  0  )  which  are  governed  by  the  level  set  equation 
(Equation  3.8).  The  steady  state  form  of  the  level  set  equation  is  the  Eikonal  equation 
(Equation  3.7)  which  says  that  the  gradient  of  the  arrival  time  surface  is  inversely 
proportional  to  the  speed  of  the  front.  For  the  two  dimensional  case,  the  stationary 
boundary  value  problem  is  given  by: 

\VT\F[x,y)  =  l  s.t.  T  =  {(a;,|/)|T(a;,|/)  =  0}  (3.9) 
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where  T  is  the  starting  position  of  the  interface.  The  first  order  hnite  difference 
discretization  form  of  the  Eikonal  equation  (Sethian,  1999b)  at  the  grid  point  (i,j)  is 
given  by: 


[max{D^^T,  0)^  +  min{D±^T,  0)^  +  maxiD^^T,  0)^  +  min{Dt^T,  0) ^ 

^  ii 


[max{max{Dj^j^T,  0),  —min{Df^T ,  0))^  + 

\ 

max{max{D~j^T,{i), —min{Df^T,{i)Y]  =  (3.10) 

Equation  3.10  is  essentially  a  quadratic  equation  for  the  value  at  each  grid  point 
(assuming  that  values  at  the  neighboring  nodes  are  known).  An  iterative  algorithm  for 
computing  the  solution  to  Equation  3.10  was  introduced  by  Ruoy  and  Tourin  (1992). 
EMM  is  based  on  the  observation  that  the  upwind  difference  structure  of  Equation 
3.10  means  that  the  information  propagates  “one  way”,  i.e.  from  the  smaller  values 
of  T  to  the  larger  values.  Therefore,  EMM  rests  on  solving  Equation  3.10  by  building 
the  solution  outward  from  the  smallest  time  value  T.  The  front  is  swept  ahead  in  an 
upwind  manner  by  considering  a  set  of  points  in  a  narrow  band  around  the  existing 
front  and  bringing  new  points  into  the  narrow  band  structure.  The  fast  marching 
algorithm  is: 

1.  Initialize 

(a)  Alive  points:  Let  A  be  the  set  of  all  grid  points  (i,j)  on  the  starting  position 
of  the  interface  T;  set  Ty  =  0  for  all  points  in  A. 

(b)  Narrow  Band  points:  Let  the  Narrow  Band  be  the  set  of  all  grid  points 
(i,j)  in  the  immediate  neighborhood  of  A;  set  Tij  =  for  all  points  in  the 
Narrow  Band  where,  d  is  the  grid  separation  distance. 

(c)  Far  Away  points:  Let  the  Far  Away  region  be  the  set  of  all  remaining  grid 
points  (i,j);  set  T^j  =  oo  for  all  points  in  the  Far  Away  region. 
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2.  Marching  Forward 


(a)  Begin  Loop:  Let  {iminjmin)  be  the  point  in  the  Narrow  Band  with  the 
smallest  value  for  T. 

(b)  Add  the  point  {imindmin)  to  A;  remove  it  from  the  Narrow  Band. 

(c)  Tag  as  neighbors  any  points  {imin  -  Ijmm),  {imin  +  IJmin),  {imindmin  -  1), 
dmindmin  +  1)  that  are  either  in  the  Narrow  Band  or  the  Far  Away  region. 
If  the  neighbor  is  in  the  Far  Away  region,  remove  it  from  that  list  and  add 
it  to  the  Narrow  Band. 

(d)  Recompute  values  of  T  at  all  neighbors  in  accordance  with  Equation  3.10. 
Select  the  largest  possible  solution  to  the  quadratic  equation. 

(e)  Return  to  the  top. 

Here  are  some  properties  of  the  fast  marching  algorithm.  The  smallest  value  in  the 
Narrow  Band  is  always  correct.  Other  Narrow  Band  or  Far  Away  points  with  larger 
values  of  T  cannot  affect  the  smallest  value.  Also,  the  process  of  recomputing  T  values 
at  the  neighboring  points  cannot  give  a  value  smaller  than  any  of  the  accepted  value 
at  Alive  points,  since  the  correct  solution  is  obtained  by  selecting  the  largest  possible 
solution  to  the  quadratic  equation  (Equation  3.10).  Thus  the  algorithm  marches 
forward  by  selecting  the  minimal  T  value  in  the  Narrow  Band  and  recomputing  the 
values  of  T  at  all  neighbors  in  accordance  with  Equation  3.10. 

The  key  to  an  efficient  version  of  the  algorithm  lies  in  Ending  a  fast  way  to  lo¬ 
cate  the  grid  point  in  the  Narrow  Band  with  the  minimum  value  for  T.  To  do  so, 
the  heapsort  algorithm  (Williams,  1964,  Sedgewick,  1988)  with  backpointers  is  often 
implemented  and  it  is  the  algorithm  we  used  here.  This  sorting  algorithm  generates 
a  “complete  binary  tree”  with  the  property  that  the  value  at  any  given  parent  node 
is  less  than  or  equal  to  the  value  at  its  child  node.  Heap  is  represented  sequentially 
by  storing  a  parent  node  at  the  location  k  and  its  child  at  locations  2k  and  2A;  -|-  1. 
The  member  having  the  smallest  value  is  stored  at  the  location  k  =  1. 

All  Narrow  Band  points  are  initially  sorted  in  a  heapsort.  The  fast  marching  al¬ 
gorithm  works  by  first  finding,  and  then  removing,  the  member  corresponding  to  the 
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smallest  T  value  from  the  Narrow  Band  which  is  followed  by  one  sweep  of  DownHeap 
to  ensure  that  the  remaining  elements  satisfy  the  heap  property.  The  DownHeap 
operation  moves  the  element  downwards  in  the  heap  till  the  new  heap  satishes  the 
heap  properties.  Far  Away  neighbors  are  added  to  the  heap  using  the  Insert  operation 
which  increases  the  heap  size  by  one  and  brings  the  new  element  to  its  correct  heap 
location  using  the  UpHeap  operation.  The  UpHeap  operation  moves  the  element  up¬ 
wards  in  the  heap  till  the  new  heap  satishes  the  heap  properties.  The  updated  values 
at  the  neighbor  points  obtained  from  Equation  3.10  are  also  brought  to  the  correct 
heap  location  by  performing  the  UpHeap  operation. 

Numerics  and  operation  count  for  the  FMM:  MATLAB  code  has  been  de¬ 
veloped  for  Objective  Analysis  using  the  Fast  Marching  Method.  Once  again,  for 
estimating  the  optimal  distance,  the  scalar  speed  function  F  is  set  to  0  for  the  grid 
points  on  land  and  1  for  the  grid  points  on  the  water.  FMM  has  a  signihcantly  lower 
operation  count  of  0(N^  Log  N)  for  grid  points  (Sethian,  1999b).  Thus,  it  is  a 
computationally  inexpensive  technique  as  compared  to  the  Level  Set  Method. 

The  Fast  Marching  Method,  as  discussed  above,  is  an  efficient  way  to  obtain  the 
correlation  between  two  locations  by  selecting  the  optimal  path.  The  length  of  the 
optimal  path  computed  using  FMM  or  LSM  can  then  be  used  for  setting  up  the 
covariance  matrix  using  the  analytical  correlation  function  (Equation  2.4). 


3.2  Objective  Analysis  based  on  using  Stochasti¬ 
cally  Forced  Differential  Equations  (SDE’s) 

The  use  of  Euclidean  distance  in  the  held  covariance  computed  from  the  isotropic 
correlation  function  is  not  applicable  in  coastal  regions  since  the  complex  coastline 
constraints,  e.g.  there  should  be  no  direct  relationship  across  landforms  (islands, 
peninsulas  etc.),  need  to  be  accounted  for.  The  approach  discussed  in  this  section 
represents  the  held  and  its  coastline  constraints  by  a  partial  diherential  equation 
subject  to  stochastic  forcing.  The  central  idea  of  this  approach,  which  is  based  on 
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using  stochastically  forced  differential  equations  (SDE),  is  the  numerical  construction 
of  a  held  covariance  such  that  it  accounts  for  the  coastal  constraints.  The  underlying 
held  variability  is  represented  as  an  outcome  of  a  stochastic  process  using  a  SDE  and 
the  stochasticity  represents  the  uncertainty  in  this  diherential  equation.  For  example, 
the  stochastically  forced  Helmholtz  equations  in  1-D  and  2-D  in  space  for  the  held  ip 
in  an  unbounded  domain  (Balgovind  et  ah,  1983)  are  associated  with  the  following 
covariance  functions  respectively: 


d^iP 


-  k'^ijj  =  e(x)  C^^(r)  =  (1  +  kr)i 


dx^ 

V^ip  -  k'^ip  =  e{x,  y)  C^^{r) 


krKi{kr)  ~  +  8^)  ^ 


oo 


(3.11) 


where,  Ki  is  the  Bessel  function  of  the  second  kind.  The  process  noise  e  is  a  random 
disturbance  with  mean  0,  standard  deviation  1  and  has  no  spatial  correlation.  Also, 
the  length  scale  corresponds  to  the  inverse  of  the  SDE  parameter  (k).  By  the  way, 
Denman  and  Freeland  (1985)  have  proposed  other  correlation  functions  which  can 
also  be  linked  to  the  appropriate  SDE’s. 


A  major  advantage  of  this  SDE  approach  is  that  the  held-to-held  covariance 
Cor(x,x)  can  be  computed  numerically  from  the  discretized  SDE  along  with  ap¬ 
propriate  boundary  conditions  (i.e.  no  hux  boundary  condition  across  islands)  to 
directly  account  for  the  coastline  constraints  (Lynch  and  McGillicuddy,  2001).  The 
discretization  of  SDE  equations  (Equation  3.11)  or  any  other  differential  operator  on 
a  hnite  element  grid  leads  to  the  matrix  form: 

l^l»}  =  {e}  (3.12) 

All  coastline  constraints  are  then  incorporated  automatically  in  the  discretization 
Equation  (3.12).  Since  [Gee]  =  [.f],  the  covariance  matrices  for  held-to-held  points 
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and  field-to-data  points  are  directly  obtained  from  Eqnation  3.12  and  given  by: 

Cor(x,x)  =  [A]-^[a,][A]-^  =  ([A]^[A])-^ 

Cor(x,X)  =  [A]-^[C^M]-^[Hf  =  {[Af[A])-^[Hf  (3.13) 

The  covariance  matrix  (3.13)  obtained  using  the  SDE  approach  can  be  used  along 
with  Gauss-Markov  Estimation  theory  to  perform  Objective  Analyses  in  coastal  re¬ 
gions.  A  limitation  of  this  approach  is  that  the  resulting  helds  can  be  affected  by  the 
discretization  error  associated  with  the  discretized  form  of  the  SDE.  In  fact,  we  found 
that  we  often  need  to  postprocess  (smooth  out)  the  SDE-gridded  helds  to  remove 
spurious  held  gradients.  Such  gradients,  even  when  small,  can  lead  to  spurious  ve¬ 
locities  by  aggregate  integration  in  the  vertical  for  the  estimation  of  absolute  velocity 
under  geostrophic  balance.  It  has  also  been  verihed  that  that  the  SDE  approach  is 
computationally  expensive  when  compared  to  our  new  FMM-based  methodology. 

A  similar  variant  of  the  above  methodology  represents  the  covariance  in¬ 

stead  of  the  held  {ip),  by  a  SDE  like  Helmholtz  equation  (Logoutov,  personal  commu¬ 
nication).  Spatial  variation  in  the  resulting  OA  held  are  found  to  be  more  prominent 
with  this  new  scheme.  An  heuristic  reason  is  that  this  new  representation  corresponds 
to  carrying  out  “smoothing”  using  the  Helmholtz  equation  only  once  as  compared  to 
twice  in  the  original  representation.  Both  of  these  methods,  the  SDE  specihed  for 
the  held  {ip)  and  the  SDE  specihed  for  the  covariance  (G^^)  were  implemented  in 
MATLAB,  using  guidance  from  Logoutov  (personal  communication). 

Even  though  many  diherent  SDE’s  could  be  utilized  for  mapping  a  held,  in  the 
example  that  follows,  we  selected  the  stochastically  forced  Helmholtz  equation.  First, 
the  dynamics  of  the  atmosphere  can  be  approximately  governed  on  the  time  scale  of 
a  few  days  by  a  Helmholtz-like  equation,  which  is  the  equation  for  the  conservation  of 
potential  vorticity  under  the  assumptions  of  a  quasi-geostrophic,  frictionless,  shallow 
water  model  without  topography  (Balgovind  et  ah,  1983,  Pedlosky,  1987).  Second, 
the  Helmholtz  equation  can  also  be  reduced  from  the  dihusion  or  wave  equations. 
In  these  linear  PDF’s,  if  the  solution  is  assumed  separable  in  time  and  space,  one 
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obtains  for  the  time  variation  an  ordinary  differential  equation  of  the  hrst  order. 
For  the  spatial  variations,  one  always  obtains  the  Helmholtz  equation  (Selvadurai, 
2000),  which  is  the  equation  that  would  be  used  for  spatial  mapping.  Thirdly,  the 
Helmholtz  equation  is  also  equivalent  to  the  steady  state  diffusion-reaction  equation. 
The  Helmholtz  equation  can  also  be  obtained  by  discretizing  the  diffusion  equation 
in  a  single  time  step. 

In  our  examples  in  Equation  3.11,  the  SDE  parameter  (k)  is  chosen  such  that  the 
correlation  function  corresponding  to  the  stochastically  forced  Helmholtz  equation 
best  hts  the  analytical  correlation  function  used  by  the  standard  OA  scheme  and  the 
LSM  or  FMM-based  schemes  (Section  3.1).  These  methods  are  compared  to  each 
other  and  to  the  EMM  and  LSM  schemes  in  Chapter  4  using  the  World  Ocean  Atlas, 
2005  data  in  the  sub-domain  of  Philippines  Archipelago. 


3.3  Estimation  of  the  absolute  velocity  under 

geostrophic  balance  by  minimizing  the  inter¬ 
island  transport  using  FMM 


For  ocean  flows,  which  evolve  over  long  spatial-time  scales  and  away  from  the  im¬ 
mediate  vicinity  of  the  sea-surface,  the  dominant  terms  in  the  horizontal  momentum 
equations  are  the  terms  corresponding  to  the  Coriolis  force  and  the  pressure  gradient. 
Such  a  flow  held,  where  a  balance  is  struck  between  the  Coriolis  and  the  pressure 
forces,  is  called  geostrophic.  The  thermal  wind  equations,  which  have  also  been  cen¬ 
tral  to  physical  oceanography  for  over  100  years,  are  obtained  for  geostrophic  how  by 
assuming  that  the  vertical  momentum  equation  is  approximately  given  by  hydrostatic 
balance.  The  thermal  wind  equations  are: 


j,d{pv)  dp  ,  ,d{pu)  dp 

—  =  9^ 

dz  dx  dz  dy 


(3.14) 
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where,  p  is  the  density,  u  and  v  are  the  horizontal  fluid  velocity  in  the  zonal  (x)  and 
meridional  (y)  directions  respectively,  and  /  =  2VL  sin(j)  is  the  Coriolis  parameter  for 
the  spherical  earth  rotating  at  a  rate  of  hi  at  latitude  0.  The  thermal  wind  equations 
(Equation  3.14)  when  integrated  in  the  vertical  give: 


pv{x,y,z,t) 


^dz  +  pMo(3.15) 
dy 


where,  Zq  or  the  level  of  no  motion  for  uq,  Mq  =  0  or  a  level  of  reference  for  uq;  Uq  ^  0. 

Flow  estimation  based  on  thermal  wind  balance  (Equation  3.15),  is  a  classical 
problem  in  physical  oceanography  (Wunsch,  1996).  Historically,  the  only  significant 
routine  measurements  possible  were  the  temperature,  T,  and  salinity,  S,  of  the  water 
at  various  depths.  The  equation  of  state  for  seawater  then  permits  the  estimation  of 
density  at  a  given  pressure  from  the  temperature  and  salinity  measurements.  Thus  the 
geostrophic  flow  can  be  computed  using  the  above  method  (Equation  3.15)  from  the 
shipboard  measurements  of  T  and  S  alone.  The  formulation  has  been  well  defined 
for  the  open  oceans  without  any  landforms.  For  complex  coastal  regions  having 
landforms  such  as  islands  and  peninsulas,  estimation  of  the  inter-island  transport  is 
first  required  before  proceeding  with  the  geostrophic  formulation  discussed  above. 

The  optimization  methodology  for  estimating  the  inter-island  transport,  proposed 
by  Haley  and  Lermusiaux  (2009)  is  utilized  and  discussed  below.  The  objective  of  this 
methodology  is  to  And  a  set  of  constant  values  for  the  transport  streamfunction  (T) 
along  the  island  coastlines  that  produce  a  suitably  smooth  initialization  velocity  held, 
e.g.  with  the  fewest  large  velocity  hot-spots,  i.e.  minimize  the  maximum  absolute 
velocity  in  the  initialized  geostrophic  flow  field.  The  working  assumptions  for  Haley’s 
methodology  are  listed  below: 

1.  Coastlines  in  the  given  domain  can  be  divided  into  two  distinct  subsets: 

(a) .  Set  A:  N  coastlines  along  which  the  transport  streamfunction  is  unknown, 
^  0. 

(b) .  Set  B:  M  coastlines  along  which  the  transport  streamfunction  is  known. 

2.  The  solution  for  the  transport  streamfunction  Tq  exists  for  the  case  which  includes 
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coasts  in  set  B,  but  coasts  in  set  A,  along  with  the  corresponding  interiors,  are  replaced 
by  open  ocean  (e.g.  island  sunk  to  lOm  depth). 

3.  The  difference  between  the  initial  solution  Tq  and  the  hnal  solution  T  is  not 
extremely  large.  Otherwise,  the  information  from  Tq  would  not  be  accurate  enough. 

To  contains  useful  information  like  the  relative  position  of  major  currents  to  var¬ 
ious  coastlines  and  the  effects  of  topography  on  the  flow.  Thus,  the  information  in 
To  can  be  utilized  to  estimate  the  constant  value  of  the  transport  streamfunction 
along  the  island  coastlines  by  constructing  an  optimization  functional  for  minimiz¬ 
ing  the  inter-island  transport  subject  to  weak  constraints.  Haley’s  methodology  for 
constructing  the  optimization  functional  is  now  discussed. 

The  problem  is  divided  into  three  parts  to  construct  the  optimization  functional. 
The  optimization  functional  (E)  in  the  general  form,  which  is  a  summation  of  three 
terms,  is  given  by: 


E  —  El  -|-  E2  -|-  E/3  (3.16) 

where,  Ei  is  the  minimizing  target  for  the  transport  between  all  pairs  of  the  unknown 
(Set  A)  coasts,  E2  is  the  minimizing  target  for  the  transport  between  all  pairs  of 
unknown  (Set  A)  and  known  (Set  B)  coasts  and  E^  is  the  minimizing  target  for  the 
transport  between  all  pairs  of  the  unknown  (Set  A)  coasts  and  the  open  boundaries 
of  the  domain.  These  three  terms  in  Equation  3.16  are: 

1.  Constructing  the  optimization  functional  for  minimizing  the  transport  between  all 

pairs  of  island  coastlines  with  unknown  (Set  A)  transport  streamfunction:  Let  and 

Cm  be  two  of  the  coasts  (coast  n  and  coast  m)  in  Set  A.  Tq  is  not  constrained  to 

be  a  constant  along  these  coasts.  Find  the  grid  point  on  the  coastline  n  and  the 

grid  point  on  the  coastline  m  such  that  =  argmax  |Ton(f)  —  Tom(j)|  and 

hi] 

=  Ton(f°)  —  Tom(j°)-  Here,  we  denote  Tq  at  point  b  on  coastline  a  by  Toa(&)- 
The  optimization  functional  for  minimizing  the  inter-island  transport  between 
islands  n  and  m  is  given  by  (T^^  —  T^^  —  (jT^^)^,  where,  T^^  is  the  unknown 
optimized,  constant  value  of  the  transport  streamfunction  along  the  coast  a.  Haley 
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proposes  to  weight  this  optimization  function  by  Wnm  =  where,  dnm  is  the 

minimum  distance  between  Cn  and  Cm-  However,  since  the  objective  is  to  smooth  the 
resulting  initialization  velocity  flow  held,  the  above  weighting  will  be  appropriate  if  the 
ocean  depth  is  uniform  in  between  all  pairs  of  islands.  An  alternative  weighting  along 
with  its  computational  methodology  for  non-uniform  ocean  depths  will  be  proposed 
later  in  this  section. 

2.  Constructing  the  optimization  functional  for  minimizing  the  transport  between 
all  pairs  of  island  coastlines  with  known  (Set  B)  and  unknown  (Set  A)  transport 
streamf unction:  Let  be  one  of  the  coasts  along  which  the  transport  streamfunction 

is  known  (Set  B)  and  be  the  coast  in  Set  A.  \ho  is  not  constrained  to  be 

a  constant  along  C„.  Find  the  grid  point  i'  on  the  coastline  n  such  that  [i']  = 

argmax  |'Fo„(f)  -  ^c'J  and  =  d'on(*')  -  'he;- 
[*] 

The  optimization  functional  for  minimizing  the  inter-island  transport  between 
islands  n  and  m  is  given  by  {'^c„  —  ^c'^.  ~  —  ^on(*0)^-  before, 

Haley  proposes  to  weight  these  optimization  function  by  where,  d'ff.  is 

the  minimum  distance  between  Cn  and  C'^. 

3.  Constructing  the  optimization  functional  for  minimizing  the  transport  between  all 
pairs  of  island  coastlines  with  unknown  (Set  A)  transport  streamfunction  and  the  open 
boundaries  of  the  domain:  Let  {Cl)  be  the  open  boundary,  {b}  be  the  set  of  open 
boundary  points  and  Cn  be  the  coast  in  Set  A.  Tq  is  not  constrained  to  be  a  constant 
along  Cn-  Find  the  grid  point  i"  on  the  coastline  n  and  the  grid  point  b"  on  the  open 
boundary  such  that  [i” -,h'']  =  argmax  |Ton(f)  — tl'(7"(&)|  and 

[i,b] 

Here,  we  denote  Tq  at  the  point  b  on  the  open  boundary  by 

The  optimization  functional  for  minimizing  the  inter-island  transport  between 
the  island  n  and  the  open  boundary  is  given  by  — 

before,  Haley  proposes  to  weight  this  optimization  function  by  w'ff,  = 
1/dnl  where,  d'^  is  the  minimum  distance  between  Cn  and  the  open  boundary  C'(. 

The  weighted  average  of  the  optimization  functionals  constructed  from  the  above 
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three  parts  is  given  by: 


N  N 

n=l  m=l,m^n 
M 

Y.  Kk('ic,  -  4'0„(*'))''  +  yU'tc.  -  4'0„(*"))"]  (3.17) 

fc=l 

The  minimum  of  E  is  computed  by  solving  the  standard  least  square  problem  i.e  by 
setting  gradients  with  respect  to  equal  to  zero.  Therefore,  the  solution  to  the 

optimization  problem  in  Equation  3.17  is  given  by: 

N  M  N 

TcJ  ^  2Wjm  +  ^'jk  +  ^”h]  -  Y.  = 

m=l^m^j  k=l  m=l,m^j 

N  M 

+  Y  (3-18) 

m=l,mj^j  k=l 

Equation  3.18  represents  a  system  of  N  equations  which  can  be  solved  to  obtain 
the  transport  streamfunctions  ('hcj)  along  coastlines  in  set  A.  These  streamfunction 
values,  which  smooth  the  velocity  held,  will  be  used  as  Dirichlet  boundary  conditions 
while  solving  the  geostrophic  how  equations  using  the  Temperature  and  Salinity  OA 
maps.  The  illustration  of  this  methodology  in  the  complex  domain  of  Philippines 
Archipelago  is  discussed  in  Section  4.3. 

Here,  we  discuss  an  extension  of  Equation  3.18  to  involve  new  more  suitable 
weights.  Consider  the  stream  function  (T)  for  a  two-dimensional  horizontal  how.  It 
is  dehned  such  that  the  how  velocity  can  be  expressed  as: 


u  =  {u,  v) 


X  ^  u 


1  9T  _  1  9T 

H  dy’’^  H  dx 


(3.19) 


Here,  H  is  the  ocean  depth.  The  transport  between  a  pair  of  islands  having  stream- 
function  ^jJl  and  1^2  is  given  by: 


t/’2  —  t/’i  =  /  u-hdA 

JA 


(3.20) 


where,  A  is  the  vertical  area  between  the  two  islands  and  h  is  the  unit  vector  normal 
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to  the  vertical  area.  Equation  3.19  and  3.20  suggests  that  the  appropriate  weight 
function  to  smooth  the  velocity  held  should  be  Wnm  =  iMnm)  where,  Anm  is  the 
minimum  vertical  area  along  any  path  between  the  two  islands.  The  weight  function 
proposed  by  Haley  (wnm  =  will  be  appropriate  when  the  ocean  depth  is 

uniform  in  between  all  pairs  of  islands.  Since  the  ocean  depth  is  not  uniform,  a  new 
methodology  is  required  to  compute  the  minimum  area  along  any  path  between  a 
pair  of  islands.  Using  the  Fast  Marching  Method  (EMM),  which  was  described  in 
Section  3.1.2,  is  a  very  convenient  and  efficient  way  to  compute  Anm-  Simulations 
have  been  performed  with  several  other  weight  functions  to  conhrm  that  the  proposed 
weight  function  based  on  the  minimum  vertical  area  {Anm)  is  the  most  appropriate 
for  smoothing  the  velocity  flow  held  (see  Chapter  4  for  illustration). 
Computational  details  for  obtaining  the  transport  streamfunction:  Fortran- 
90  code  from  MSEAS  (Haley,  personal  communication)  has  been  modihed  to  utilize 
the  weight  function  based  on  the  minimum  vertical  area  between  islands,  which  was 
computed  using  the  EMM,  instead  of  the  weight  function  based  on  the  minimum 
distance  between  islands.  For  obtaining  the  minimum  vertical  area,  the  scalar  speed 
function  in  the  Eikonal  Equation  (Equation  3.7)  is  chosen  to  be  F(x,y)  =  1/H(x,y). 
The  cost  of  computing  the  minimum  distance  is  equal  to  the  cost  of  computing  the 
minimum  vertical  area  using  the  EMM.  Thus,  the  computational  cost  is  independent 
of  the  choice  of  the  weight  function. 

New  methodologies  for  Objective  Analysis  in  the  multiply-connected  coastal  re¬ 
gions  were  proposed  in  this  Chapter.  The  applications  of  these  methodologies  are 
discussed  in  Chapter  4.  The  Level  Set  Method  and  the  Fast  Marching  Method  have 
been  used  for  OA  in  the  complex  domains  of  the  Philippines  Archipelago  and  Dabob 
Bay.  These  methodologies  have  been  compared  with  the  approach  proposed  by  Lynch 
and  McGillicuddy  (2001)  and  the  standard  OA.  Estimation  of  absolute  velocity  under 
geostrophic  balance  in  Philippines  Archipelago  is  also  illustrated  in  Chapter  4.  This 
is  followed  by  a  computational  study  of  properties  of  the  new  mapping  schemes  in 
Chapter  5. 
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Chapter  4 


Applications  illustrating  the  novel 
OA  methodologies 


New  methodologies  for  Objective  Analysis  in  complex  multiply-connected  coastal 
regions  were  described  in  Chapter  3.  These  new  methodologies  are  based  on  com¬ 
puting  the  optimal  path  lengths  using  the  Level  Set  Method  and  the  Fast  Marching 
Method.  These  methods  efficiently  incorporate  coastline  constraints  (e.g.  there  is  no 
direct  relationship  across  landforms).  The  above  methodologies  are  utilized  to  map 
the  temperature,  salinity  and  biological  (chlorophyll)  fields  using  a  2-staged  mapping 
scheme  in  subsets  of  the  following  regions:  Dabob  Bay  and  Philippines  Archipelago. 

Section  4.1  evaluates  the  use  of  our  new  OA  methodology  in  Dabob  Bay  and  shows 
that  it  is  more  effective  over  other  classic  distance  optimizing  algorithms  like  Bresen- 
ham’s  line  algorithm  (Bresenham,  1965).  Section  4.2  shows  a  comparison  of  the  dif¬ 
ferent  methodologies  introduced  in  Chapter  3  for  Objective  Analysis  in  a  subdomain 
of  Philippines  Archipelago.  The  estimation  of  absolute  velocity  under  geostrophic 
balance  by  minimizing  the  inter-island  transport  is  illustrated  in  Section  4.3.  Section 
4.4  discusses  approaches  for  including  non-homogeneous  dynamical  effects  in  our  new 
Objective  Analysis  schemes. 
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4.1  Objective  Analysis  in  Dabob  Bay 


Dabob  Bay  data  are  used  to  illustrate  the  effectiveness  of  the  Fast  Marching 
Method  over  other  distance  optimizing  algorithms  like  Bresenham  line  algorithm  (Bre- 
senham,  1965).  Maps  for  the  temperature  and  salinity  fields  in  a  subdomain  of  Dabob 
Bay  corresponding  to  the  spatially  irregular  data  in  Figure  C-4  are  obtained  using 
the  a.  Bresenham  line  algorithm,  b.  Averaged  Bresenham  line  algorithm,  and  c.  Fast 
Marching  Method.  The  limitation  of  Bresenham  line  algorithm,  which  is  explained  in 
Chapter  3,  is  that  the  optimal  distance  computed  using  this  method  is  discontinuous. 
This  results  in  discontinuities  in  the  covariance  and  also  in  the  resultant  field  maps. 

Figure  C-5  shows  the  temperature  and  salinity  field  maps  in  Dabob  bay  obtained 
using  large  length  scales  (Lq  =  60,  L^.  =  30)ls,  most  energetic  length  scales  (Lq  =  30, 
Lf.  =  15)me  and  observational  error  {R  =  0.25/).  Temperature  and  salinity  have 
higher  magnitudes  in  the  northern  part  of  the  western  arm  of  Dabob  bay.  The  eastern 
arm  of  Dabob  bay  has  relatively  low  temperature  and  salinity.  Effects  due  to  the 
discontinuity  in  distance  obtained  from  Bresenham  line  algorithm  is  clearly  evident 
in  Figure  C-5(top).  Numerical  fronts  having  high  temperature  and  salinity  gradients 
exist  at  the  intersection  of  the  two  arms.  Such  fronts  lead  to  numerical  problems 
in  dynamical  simulations.  The  geostrophic  velocity  obtained  using  these  field  maps 
will  be  unrealistic  and  will  have  high  magnitudes  along  these  fronts.  A  possible 
remedy,  which  reduces  the  discontinuity  effects,  is  to  smooth  the  distance  by  averaging 
distances  of  neighboring  points  (Lermusiaux  and  Haley,  personal  communication). 
The  above  averaging  technique  becomes  numerically  very  expensive  as  shown  by 
Lermusiaux  and  Haley. 

The  field  maps  obtained  using  the  averaged  Bresenham  algorithm  (Figure  C- 
5 (middle))  clearly  show  that  the  intensity  of  the  erroneous  fronts  are  reduced,  but 
they  still  exist.  Finally,  the  Fast  Marching  Method  is  used  to  compute  distances 
and  the  Objective  Analysis  field  maps  obtained  using  FMM  are  clearly  devoid  of  any 
numerical  fronts  (Figure  C-5  (bottom)).  Along  with  that,  FMM  accurately  satisfies 
coastline  constraints  and  it  is  computationally  inexpensive  compared  to  using  the 
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averaged  Bresenham  line  algorithm.  Thus,  the  use  of  FMM  is  recommended  over  the 
Bresenham  line  algorithm. 


4.2  Objective  Analysis  in  the  Philippines 
Archipelago 

This  research  study  is  motivated  by  the  Philippines  Straits  Dynamics  Experiment 
(PhilEx)  sponsored  by  the  Office  of  Naval  Research.  Novel  OA  techniques  for  such 
complex  coastal  regions  are  an  important  requirement  to  map  very  irregular  datasets 
and  initialize  simulations.  A  comparison  of  the  different  OA  methodologies  will  be 
illustrated  in  this  region.  We  compare  our  new  methods  a.  Level  Set  Method,  b. 
Fast  Marching  Method,  to  the  existing  schemes,  the  a.  Standard  OA  Method  which 
ignores  islands  and  uses  the  direct  Euclidean  distance,  b.  Stochastically  forced  Dif¬ 
ferential  Equation  approach  (SDE  specihed  for  the  held)  and  c.  Stochastically  forced 
Differential  Equation  approach  (SDE  specihed  for  the  covariance). 

A  comparison  of  these  methods  using  the  World  Ocean  Atlas,  2005  (Locarnini 
et  ah,  2006,  Antonov  et  ah,  2006)  data  for  the  temperature  and  salinity  held  maps  is 
discussed  in  Section  4.2.1.  WOA-05  data  are  data  mapped  using  ‘Levitus  climatology’ 
scheme  (see  Section  2.2)  and  is  regularly  spaced.  Regularly  spaced  WOA-05  data  is 
used  here  primarily  to  illustrate  and  discuss  the  comparison  of  the  diherent  method¬ 
ologies.  Subsequently,  synoptic  in  situ  data  is  used.  These  real  direct  ocean  data  are 
the  spatially  irregular  temperature  and  salinity  data.  Results  are  presented  in  Section 

4.2.2  (using  Melville  exploratory  cruise  data.  Global  Temperature-Salinity  Prohle  Pro¬ 
gram  -  GTSPP  (Genter,  2006)  data  and  HB2  Glimatology  for  June-July’07),  section 

4.2.3  (using  Melville  exploratory  cruise,  sgl22  and  sgl26  glider  data  for  June-July’07) 
and  section  4.2.4  (using  joint  Melville  cruise  data  for  Nov’07-Jan’08).  Results  of  the 
new  OA  methodology  based  on  FMM  is  illustrated  in  section  4.2.5  for  the  biologi¬ 
cal  held  (chlorophyll)  using  the  exploratory  cruise  data.  These  biological  OA  held 
maps  obtained  using  our  FMM  scheme  can  be  utilized  in  the  initialization  for  coupled 
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physics-biology  modeling  in  MSEAS  (Burton,  2009). 


4.2.1  Objective  Analysis  using  WOA-05  data:  Comparison 
of  the  different  OA  methodologies 

Two-dimensional  horizontal  OA  maps  for  temperature  and  salinity  fields  in  a 
subdomain  of  the  Philippines  Archipelago  corresponding  to  the  data  in  Figure  C-6 
computed  using  methodologies  proposed  in  Chapter  3  are  shown  in  Figures  C-7  - 
C-16.  Figures  C-7  (LSM),  C-8  (FMM),  C-9  (Standard),  C-10  (SDE  specihed  for  the 
covariance)  and  C-11  (SDE  specihed  for  the  held)  show  the  temperature  held  maps, 
while  Figures  C-12  (LSM),  C-13  (FMM),  C-14  (Standard),  C-15  (SDE  specihed  for 
the  covariance)  and  C-16  (SDE  specihed  for  the  held)  show  the  salinity  held  maps. 
Depths  shown  are  Om,  40m,  200m,  450m,  1000m  and  3000m.  Large  length  scales 
{Lq  =  540,  Le  =  180)ls  and  most  energetic  length  scales  (Lq  =  180,  Lg  =  60)me  are 
used  with  an  observational  error  covariance  R  =  0.25J.  For  the  SDE  approach,  the 
SDE  parameter  k  =  1/200  and  the  observational  error  [R  =  0.25/)  are  used. 

The  OA  held  maps  from  all  methods  indicate  that  the  Philippines  Sea  and  the 
region  near  Palawan  island  is  warmer  than  the  rest  of  the  region  at  the  surface  (Om). 
The  region  south  of  the  Sulu  sea  around  the  Sulu  Archipelago  has  relatively  lower 
temperature.  At  levels  below  500m,  there  is  a  signihcant  diherence  in  the  temperature 
of  the  Sulu  sea  (warm)  as  compared  to  the  rest  of  the  region  (cold)  (Gamo  et  ah, 
2007,  Gordon,  2009).  These  temperature  helds  clearly  show  that  direct  correlation 
across  landforms  are  weak.  Similar  observations  can  be  made  for  Salinity.  Salinity  in 
the  Sulu  Sea  and  South  China  Sea  is  lower  than  the  salinity  in  the  rest  of  the  region 
at  the  surface  (Om).  At  levels  below  500m,  the  salinity  in  the  Sulu  sea  is  signihcantly 
lower  as  compared  to  the  rest  of  the  region.  These  salinity  helds  further  support  the 
hypothesis  that  direct  correlation  across  landforms  are  weak. 

The  comparison  of  the  temperature  held  maps  and  the  salinity  held  maps  ob¬ 
tained  using  diherent  methods  at  level  1000m  is  shown  in  Figures  C-17  and  C-18, 
respectively.  The  methods  based  on  LSM,  FMM  and  SDE  clearly  satisfy  coastline 
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constraints.  The  data  in  the  Sulu  Sea,  which  has  high  temperature  and  low  salinity 
compared  to  the  remaining  region,  does  not  have  any  influence  on  the  held  outside 
the  Sulu  Sea  since  the  two  regions  are  not  connected  by  water.  On  the  other  hand, 
the  standard  OA  does  not  satisfy  coastline  constraints.  Thus  the  data  outside  the 
Sulu  Sea,  where  the  temperature  is  low  and  salinity  is  high,  is  correlated  to  the  held 
inside  the  Sulu  Sea.  This  is  undesirable  since  the  direct  relationship  across  landforms 
is  at  best  very  weak.  This  leads  to  spurious  high  temperature  and  salinity  gradients 
in  the  Sulu  Sea,  which  will  lead  to  problems  for  the  estimation  of  geostrophic  flow. 
Differences  between  temperature  field  maps  and  salinity  field  maps  obtained  using 
the  Fast  Marching  Method  and  using  other  OA  methods  at  level  1000m  are  shown 
in  Figures  C-19  and  C-20,  respectively.  These  plots  show  that  there  is  a  very  small 
difference  in  held  maps  obtained  using  the  FMM  and  LSM.  The  difference  is  larger 
between  field  maps  obtained  using  the  FMM  and  SDE  approach.  This  is  because  the 
analytical  correlation  function  corresponding  to  the  stochastically  forced  Helmholtz 
equation,  which  is  used  in  the  SDE  approach,  is  different  from  the  analytical  correla¬ 
tion  function  in  the  FMM.  The  difference  between  field  maps  obtained  using  the  FMM 
and  standard  OA  are  significantly  large  because  standard  OA  does  not  incorporate 
coastline  constraints. 

The  held  maps  obtained  by  LSM  (Figures  C-7,  C-12)  and  FMM  (Figures  C-8, 
C-13)  are  almost  identical,  but  the  FMM  has  a  significantly  lower  computational 
cost.  While  LSM  constructs  the  solution  by  iterating  towards  the  solution,  FMM  is 
based  on  the  direct  construction  of  the  stationary  solution  as  described  in  Section  3.1. 
There  is  a  very  small  difference  in  the  held  obtained  using  LSM  and  FMM  because 
FMM  exactly  constructs  the  solution  of  the  discretized  Eikonal  equation  whereas  LSM 
computes  the  solution  within  a  desired  tolerance  limit.  So,  FMM  is  more  accurate 
and  less  expensive  compared  to  LSM.  Thus,  our  OA  methodology  based  on  FMM 
should  clearly  be  preferred  over  our  methodology  based  on  LSM. 

The  SDE  approach  satisfies  coastline  constraints,  but  the  discretization  errors 
in  SDE  are  significant  and  this  results  in  prominent  spatial  variations  in  the  tem¬ 
perature  and  salinity  fields.  The  impact  of  such  huge  spatial  variations  on  the 
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geostrophic  flow  velocity  is  not  good,  and  often  additional  smoothing  has  to  be  em¬ 
ployed  (post-processing)  after  obtaining  the  OA  helds  using  the  SDE  approach.  Such 
post-processing  is  not  required  for  our  FMM-based  scheme.  The  SDE  approach  can 
be  implemented  by  specifying  the  SDE  for  the  held  or  by  specifying  it  for  the  co- 
variance  (Logoutov,  Lermusiaux,  personal  communication).  If  the  SDE  is  specihed 
for  the  held  (Figure  C-11,  C-16)  as  opposed  to  the  covariance  (Figure  C-10,  C-15), 
spatial  variation  in  the  held  will  be  less  prominent.  Specifying  SDE  (say  Helmholtz 
equation)  for  the  held  is  equivalent  to  carrying  out  the  smoothing  twice  by  using  the 
Helmholtz  equation.  Thus,  specifying  the  SDE  for  the  held  will  be  more  expensive 
than  specifying  the  SDE  for  the  covariance,  but  this  will  make  the  spatial  variation 
in  the  held  less  prominent  and  it  will  reduce  the  need  for  post-processing.  Finally,  we 
conhrmed  that  the  computational  time  required  by  the  SDE  approach  is  higher  than 
that  of  EMM.  Thus,  EMM  appears  to  be  the  best  among  all  the  methodologies  dis¬ 
cussed  in  Chapter  3  for  Objective  Analysis  in  coastal  regions.  Therefore,  in  Sections 
4.2.2,  4.2.3,  4.2.4  and  4.2.5,  we  show  and  discuss  results  of  our  FMM-based  Objective 
Analysis  scheme  for  mapping  spatially  irregular  data. 

4.2.2  Objective  Analysis  for  Summer  2007:  Melville 

exploratory  cruise,  GTSPP  and  HB2  Climatology  data 

The  coarse  WOA-05  data  used  in  Section  4.2.1  is  regularly  spaced  and  is  already 
mapped  using  the  ‘Levitus  Climatology’  mapping  scheme.  The  data  used  now  is 
sampled  in  situ  and  is  irregularly  spaced.  They  were  collected  from  the  Melville  Ex¬ 
ploratory  cruise,  the  Global  Temperature-Salinity  Prohle  Program  -  GTSPP  (Center, 
2006)  and  the  HB2  Climatology  for  the  June-July’07  period.  The  data  location  plot 
is  shown  in  Figure  C-21.  Large  length  scales  (Lq  =  1080,  Lg  =  360)ls,  most  energetic 
length  scales  (Lq  =  270,  =  90)me  and  observational  error  (i?  =  0.2/)  are  used. 

The  temperature  and  salinity  held  maps  obtained  using  the  FMM-based  OA 
scheme  at  the  surface  (Om)  in  Figure  C-22  clearly  show  that  coastline  constraints 
are  appropriately  incorporated,  since  the  warm  region  in  the  west  of  Luzon  island  is 
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uncorrelated  with  the  region  on  the  east  of  Luzon  island.  The  same  holds  true  for 
salinity,  since  the  data  in  the  low  salinity  region  i.e.  west  of  Luzon  island  does  not 
have  a  signihcant  impact  on  the  held  in  the  east  of  Luzon  island.  It  is  also  observed 
that  local  effects  and  wind  patterns  lead  to  higher  temperature  and  salinity  in  the 
Visayan  Sea  while  the  temperature  and  salinity  in  the  Bohol  Sea  remain  low. 

Again,  it  should  be  noted  that  the  data  from  the  HB2  climatology  is  mapped 
data.  Section  4.2.3  and  4.2.4  will  discuss  the  examples  of  Objective  Analysis  using 
the  exploratory  cruise  data  alone. 

4.2.3  Objective  Analysis  for  Summer  2007:  Melville 
exploratory  cruise,  sgl22  and  sgl26  glider  data 

The  data  used  in  this  example  is  collected  from  the  Melville  exploratory  cruise, 
sgl22  and  sgl26  gliders  for  the  June-July’07  period.  The  data  location  plot  is  shown 
in  Figure  C-23.  Since  the  data  is  available  only  in  a  small  region  of  the  Philippines 
Archipelago  near  islands.  Objective  Analysis  maps  are  computed  in  a  portion  of  the 
regular  Philex  domain.  Large  length  scales  (Lq  =  1080,  Lg  =  360) ls,  most  energetic 
length  scales  (Lq  =  270,  Lg  =  90)me  and  observational  error  {R  =  0.2/)  are  used. 
The  temperature  and  salinity  held  maps  obtained  using  the  methodology  based  on  the 
Fast  Marching  Method  are  shown  in  Figures  C-24  and  C-25,  respectively  at  depths  of 
Om,  40m,  200m,  450m,  1000m  and  3000m.  Once  again,  these  maps  clearly  indicate 
that  coastline  constraints  are  appropriately  satished.  At  depths  of  Om  and  40m,  the 
warm  region  in  the  west  of  Luzon  island  is  uncorrelated  with  the  region  on  the  east 
of  Luzon  island.  The  warm  Sibuyan  and  Visayan  Seas  can  be  distinguished  from 
the  relatively  cold  Bohol  Sea.  At  depths  of  450m  and  1000m,  the  data  in  the  warm 
Sulu  sea  and  Bohol  Sea  does  not  have  any  impact  on  the  remaining  regions,  clearly 
suggesting  that  there  is  no  direct  relationship  across  landforms.  Similar  observations 
are  made  for  the  salinity.  At  depths  of  Om  and  40m,  the  low  salinity  region  in  the 
west  of  Luzon  island  is  uncorrelated  with  the  region  on  the  east  of  Luzon  island. 

We  now  compare  helds  in  Snmmer  2007  from  Melville  exploratory  crnise,  sgl22 
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and  sgl26  glider  data  (Figures  C-24  and  C-25)  with  fields  from  from  Melville  ex¬ 
ploratory  cruise,  GTSPP  and  HB2  Climatology  data  (Figure  C-22)  (Section  4.2.3). 
This  comparison  illustrates  differences  in  fields  obtained  using  different  datasets  for 
Summer  2007.  There  is  a  good  comparison  between  the  two  fields  near  islands.  While 
the  exploratory  cruise  and  glider  data  is  available  in  a  small  region  around  islands,  the 
GTSPP  and  HB2  Climatology  data  is  also  available  away  from  islands.  Thus  away 
from  the  island,  fields  from  Melville  exploratory  cruise,  GTSPP  and  HB2  Climatology 
data  (Figure  C-22)  will  be  more  accurate. 

Section  4.2.4  will  now  discuss  the  example  of  Objective  Analysis  using  the  joint 
Melville  cruise  data  in  the  winter  season  (Nov’07-Jan’08). 

4.2.4  Objective  Analysis  for  Winter  2008:  Melville  joint 
cruise  data 

The  data  used  in  this  example  is  obtained  from  the  joint  Melville  cruise  for  the 
Nov’07-Jan’08  period.  The  data  location  plot  is  shown  in  Figure  C-26.  Once  again, 
since  data  are  available  in  a  small  region  of  the  Philippines  Archipelago  near  islands, 
maps  are  obtained  in  a  smaller  subsection  of  the  regular  Philippines  region.  Large 
length  scales  (Lq  =  1080,  =  360)ls,  most  energetic  length  scales  (Lq  =  270, 

Lf,  =  90)me  and  observational  error  (i?  =  0.2J)  are  used  for  the  OA  field  maps. 
The  temperature  and  salinity  field  maps  obtained  using  the  FMM-based  scheme  are 
shown  in  Figures  C-27  and  C-28,  respectively.  Depths  shown  are  Om,  40m,  200m, 
450m,  1000m  and  3000m.  Once  again,  at  depths  of  Om  and  40m,  the  warm  region 
in  the  west  of  Luzon  island  is  uncorrelated  with  the  region  on  the  east  of  Luzon 
island.  At  depths  of  450m  and  1000m,  the  data  in  the  warm  Bohol  Sea  does  not 
have  any  impact  on  the  remaining  regions,  clearly  suggesting  that  there  is  no  direct 
relationship  across  landforms.  Similar  observations  are  made  for  salinity.  At  depths 
of  Om  and  40m,  the  low  salinity  region  in  the  west  of  Luzon  island  is  uncorrelated 
with  the  region  in  the  east  of  Luzon  island. 

We  now  compare  fields  in  Winter  2008  from  Melville  joint  cruise  data  with  fields  in 
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Summer  20007  from  Melville  exploratory  cruise,  sgl22  and  sgl26  glider  data  (Section 
4.2.3).  It  is  clearly  evident  that  the  difference  in  temperature  during  Winter  2008 
and  Summer  2007  is  more  near  the  ocean  surface.  Beyond  the  depth  of  200m,  the 
difference  is  signihcantly  less  and  the  same  inference  is  valid  for  salinity  as  well.  At 
surface  (Om),  the  temperature  in  the  Sulu  sea  is  nearly  the  same  for  both  Summer 
2007  and  Winter  2008.  But  the  temperature  near  Luzon  island  is  signihcantly  lower 
during  Winter  2007  than  the  temperature  during  Summer  2007. 

Section  4.2.5  illustrates  the  application  of  our  FMM-based  scheme  for  biological 
held  (chlorophyll). 

4.2.5  Objective  Analysis  for  biological  field  (chlorophyll) 

Application  of  our  new  FMM-based  scheme  for  the  biological  held  (chlorophyll)  is 
illustrated  here  using  Exploratory  cruise  Summer  2007  data.  The  biological  OA  held 
map  obtained  using  FMM  can  be  utilized  in  the  initialization  for  coupled  physics- 
biology  modeling  studies  (Burton,  2009).  The  data  location  plot  is  shown  in  Figure 
C-29.  Large  length  scales  (Lq  =  1080,  L^.  =  360)ls,  most  energetic  length  scales 
(Lq  =  270,  Le  =  90)me  and  observational  error  {R  =  0.2/)  are  used  for  the  OA 
held  maps  at  depths  of  Om,  40m,  200m,  450m,  1000m,  3000m.  The  chlorophyll 
maps  computed  using  our  FMM-based  scheme  are  shown  in  Figure  C-30  at  depths 
of  Om,  10m,  50m,  100m,  150m,  200m.  The  concentration  of  biological  helds  like 
chlorophyll,  phytoplankton  and  zooplankton  is  substantial  only  near  the  surface  due 
to  the  presence  of  sunlight.  Therefore,  the  coupled  physics-biology  modeling  studies 
are  usually  carried  up  to  the  depth  of  200m. 

The  chlorophyll  concentration  is  maximum  near  islands.  Away  from  islands,  it 
approaches  the  mean  data  value.  At  depth  of  Om  and  10m,  the  maximum  chlorophyll 
concentration  is  observed  in  the  south  of  the  Visayan  sea  and  in  the  Bohol  Sea.  At  a 
depth  of  50m,  the  chlorophyll  concentration  in  the  south  of  the  Visayan  sea  and  in  the 
Bohol  Sea  remains  signihcant.  The  maximum  chlorophyll  concentration  is  observed 
in  the  north  of  Palawan  island. 
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4.3  Estimation  of  the  absolute  velocity  under 

geostrophic  balance  in  Philippines  Archipelago 

Estimation  of  absolute  velocity  under  geostrophic  balance  in  the  Philippines  Archi¬ 
pelago  is  illustrated  in  this  section.  The  algorithm  for  minimizing  the  inter-island 
transport  (Chapter  3)  is  utilized  for  computing  a  smooth  geostrophic  velocity  flow 
field. 

We  proposed  to  utilize  weight  functions  based  on  the  minimum  vertical  area  along 
each  pair  of  islands  in  the  algorithm  for  minimizing  the  inter-island  transport.  The 
estimation  of  the  minimum  vertical  area  has  been  carried  out  using  the  EMM  by 
specifying  the  scalar  speed  function  in  the  Eikonal  equation  (Equation  3.7)  as  F(x,y) 
=  1/H(x,y),  where  H  is  the  ocean  depth.  The  temperature  and  salinity  data  are 
from  the  World  Ocean  Atlas  2005  (Figure  C-6).  They  are  mapped  using  our  FMM- 
based  OA  scheme  (Figures  C-8  and  C-13)  and  the  SDE  approach  (Figures  C-11  and 
C-16),  with  the  stochastically  forced  Helmholtz  equation  employed  for  the  field.  The 
streamfunction  and  velocity  fields  (at  depths  Om,  100m,  1000m)  using  the  maps  from 
our  FMM-based  scheme  are  shown  in  Figure  C-31.  These  plots  show  a  very  good 
comparison  with  the  streamfunction  and  velocity  obtained  using  the  temperature 
and  salinity  held  maps  based  on  the  stochastically  forced  Helmholtz  equation,  which 
are  shown  in  Figure  C-32.  These  maps  suggest  that  the  velocity  is  maximum  in  the 
Mindoro  strait,  near  the  Mindanao  island  and  in  the  Balabac  strait.  At  lower  depths, 
the  velocity  remain  high  in  the  Mindoro  strait  and  near  the  Mindanao  island.  There  is 
a  large  inter-island  transport  across  the  Mindoro  strait  since  the  vertical  area  between 
the  Mindoro  and  Palawan  island  is  very  large. 

Haley’s  methodology  utilizes  weight  functions  based  on  the  minimum  inter-island 
distance  which  can  be  obtained  using  the  FMM  by  specifying  the  scalar  speed  function 
in  the  Eikonal  equation  (Equation  3.7)  as  1  for  the  sea  points  and  0  for  the  land 
points.  The  streamfunction  and  velocity  fields  (at  depths  Om,  100m,  1000m)  are 
obtained  using  the  weight  functions  based  on  the  minimum  inter-island  distance.  The 
streamfunction  and  velocity  helds  using  held  maps  from  our  FMM-based  scheme  are 
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shown  in  Figure  C-33.  Once  again,  these  plots  show  a  very  good  comparison  with  the 
streamfunction  and  velocity  obtained  using  the  temperature  and  salinity  held  maps 
based  on  the  stochastically  forced  Helmholtz  equation,  which  are  shown  in  Figure  C- 
34.  These  maps  also  suggest  that  the  velocity  is  maximum  in  the  Mindoro  strait,  near 
the  Mindanao  island  and  in  the  Balabac  strait.  But  the  velocity  estimated  in  Balabac 
strait  is  very  high  and  is  clearly  not  acceptable.  Such  high  velocity  is  obtained  due 
to  the  inaccurate  computation  of  inter-island  transport. 

These  results  clearly  show  that  the  weight  functions  based  on  the  minimum  vertical 
area  will  produce  smooth  geostrophic  flow  held  with  the  least  velocity  hot  spots. 


4.4  Non-homogeneous  dynamical  effects  in  Objec¬ 
tive  Analysis 

We  satished  coastline  constraints  using  the  FMM/LSM-based  schemes,  which 
solve  the  Eikonal  equation  (Equation  3.7)  to  obtain  the  length  of  the  optimal  path.  By 
appropriately  specifying  the  scalar  speed  function  F  in  the  Eikonal  equation,  which 
simply  says  that  the  gradient  of  the  arrival  time  function  is  inversely  proportional 
to  the  speed,  the  optimal  path  length  can  be  obtained.  The  scalar  speed  function  is 
specihed  as  1  for  the  sea  points  and  0  for  the  land  points  in  the  illustrations  in  Section 
4.1  and  4.2.  Additional  effects  due  to  the  ocean  bathymetry  and  the  dynamics  can 
also  be  incorporated  in  the  Objective  Analysis  by  appropriately  modifying  the  scalar 
speed  function  {F)  or  the  choice  of  the  length  scales. 

Such  dynamical  effects  can  be  explained  in  the  context  of  continental  a  shelf.  The 
scales  in  a  continental  shelf  are  not  uniform.  The  correlation  between  any  two  points 
in  the  ocean  will  depend  on  the  optimal  path  between  the  points  and  the  non-uniform 
length  scales  on  the  optimal  path.  While  the  optimal  path  can  be  obtained  using  the 
FMM  or  LSM,  the  utilization  of  all  the  length  scales  on  the  optimal  path  has  to 
be  appropriately  done.  It  can  be  argued  that  the  dynamics  will  be  governed  by  the 
smallest  length  scale  on  the  optimal  path,  since  the  smallest  length  scale  will  govern 
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the  correlation  between  the  two  points.  Consider  the  example  of  a  continental  shelf 
shown  in  Figure  C-35.  In  this  hgure,  the  scales  in  y-direction  are  large  in  both  the 
regions:  region  A  and  region  B.  These  regions  are  separated  by  a  front  having  small 
length  scales.  Therefore,  the  correlation  between  any  two  points  in  the  same  region 
will  be  governed  by  the  large  length  scale  and  the  correlation  between  the  points  in 
different  regions  will  be  governed  by  the  small  length  scale.  Thus,  the  effect  of  non- 
homogeneous  length  scales  can  be  incorporated  by  choosing  the  smallest  scale  on  the 
optimal  path. 

Effects  due  to  the  bathymetry  can  be  incorporated  by  specifying  the  scalar  speed 
function  at  a  given  depth  as  the  ratio  of  the  bathymetry  at  the  grid  point  and  the 
depth  at  which  the  OA  is  carried  out.  This  will  ensure  that  even  if  the  two  points 
in  the  domain  are  separated  by  land  at  a  certain  depth,  the  correlation  between 
those  points  at  that  depth  is  not  necessarily  zero.  Since  the  land  may  not  extend 
to  a  significant  height  above  the  depth  at  which  the  OA  is  carried  out,  it  will  be 
inappropriate  to  specify  the  scalar  speed  function  as  zero  for  such  a  pair  of  points. 
This  is  illustrated  in  Figure  C-36  for  the  Philippines  Archipelago  at  a  depth  of  450m. 
This  hgure  shows  the  scalar  speed  functions  for  the  general  case  and  for  the  case  in 
which  effects  due  to  the  bathymetry  are  incorporated.  The  temperature  and  salinity 
helds  obtained  for  these  different  cases  are  also  compared.  The  difference  is  clearly 
visible  in  the  Bohol  Sea.  In  the  general  case,  the  data  in  the  Sulu  Sea  is  uncorrelated 
with  the  held  in  the  Bohol  Sea.  But  the  use  of  modihed  scalar  speed  function  shows 
that  the  data  in  the  Sulu  Sea  has  an  ahect  on  the  held  in  the  Bohol  Sea.  The  scalar 
speed  function  can  also  be  modihed  using  other  non-linear  functions  which  depend 
on  the  bathymetry  at  the  grid  point  and  the  depth  at  which  the  OA  is  carried  out. 

This  discussion  of  the  modihcation  of  the  length  scales  or  the  scalar  speed  func¬ 
tion  can  similarly  be  altered  in  diherent  situations  to  include  the  dynamical  ehects 
(for  example,  conservation  of  the  potential  vorticity)  in  the  Objective  Analysis.  This 
concludes  the  demonstration  of  the  new  OA  methodologies  and  the  methodology  for 
obtaining  the  geostrophic  flow  velocities  in  complex  coastal  regions.  The  computa¬ 
tional  details  of  the  OA  methodologies  will  be  discussed  in  Chapter  5. 
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Chapter  5 


Computational  Analysis  and 
Derivations 


Computational  studies  of  properties  of  the  new  mapping  schemes  are  carried  out 
in  this  Chapter.  Section  5.1  introduces  the  sequential  processing  of  observations 
for  mapping  irregular  data  using  our  new  OA  schemes.  Sequential  processing  re¬ 
duces  computational  costs  and  it  also  allows  to  estimate  the  impact  of  the  individual 
data.  We  introduce  dehnitions  of  convex,  simply-connected  and  multiply-connected 
domains  here.  A  domain  is  said  to  be  convex  if  for  every  pair  of  points  within  the 
domain,  every  point  on  the  straight  line  segment  that  joins  them  is  also  within  the 
domain.  A  domain  is  said  to  be  simply-connected  if  any  closed  curve  within  it  can  be 
continuously  shrunk  to  a  point  without  leaving  the  domain.  A  domain  which  is  not 
simply-connected  is  called  multiply-connected.  Section  5.2  introduces  the  Wiener- 
Khinchin  and  Bochner  theorems  for  positive  dehnite  correlation  functions.  Positive 
dehnite  correlation  functions  can  generate  a  positive  dehnite  covariance  matrix  for 
a  simply-connected  convex  domain  using  the  Euclidean  distance.  Our  mapping  of 
observations  is  linked  to  the  Kalman  Filter’s  state  update.  It  is  known  that  the 
Kalman  Filter  encounters  divergence  problems  if  the  covariance  matrix  becomes  neg¬ 
ative  due  to  numerical  issues  (Brown  and  Hwang,  1997).  Some  useful  techniques  to 
counter  these  divergence  problems  are  discussed  in  Section  5.3.  The  comparison  of 
computational  costs  for  different  schemes  is  made  in  Section  5.4. 
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The  Wiener-Khinchin  and  Bochner  theorems  are  valid  for  the  background  er¬ 
ror  covariance  matrix,  Cor(x,x)  computed  using  the  Euclidean  distance  for  simply- 
connected  convex  domains.  For  complex  coastal  regions,  Cor(x,x)  may  not  necessarily 
be  positive  dehnite  due  to:  a.  numerical  error  in  the  computation  of  the  optimal  path 
length  using  FMM/LSM  b.  the  presence  of  landforms.  This  may  lead  to  divergence 
problems  for  the  held  mapping  based  on  the  FMM/LSM  scheme  in  complex  coastal 
regions.  Such  divergence  problems  are  illustrated  using  the  WOA-05  data  (Spliced 
February  and  Winter  Climatology)  shown  in  Figure  C-37.  The  held  maps  obtained 
using  our  FMM-based  scheme  (one-scale)  with  length  scales  (Lq  =  540,  Lg  =  180) 
and  length  scales  {Lq  =  1080,  Lg  =  360)  are  shown  in  Figure  C-38.  Fields  obtained 
using  the  larger  scales  clearly  show  divergence  problems  near  the  Palawan  island. 
Such  problems  are  not  encountered  when  the  smaller  length  scales  are  used.  Specif¬ 
ically,  questions  which  motivate  our  research  in  Section  5.5  and  5.6  are:  a.  What 
are  the  computational  errors  in  optimal  path  lengths  computed  using  the  FMM/LSM 
and  how  can  they  be  reduced?  b.  What  are  the  computational  issues  including  non¬ 
positive  dehnite  covariances  that  arise  in  a  multiply-connected  coastal  domain  and 
how  can  they  be  remedied?  A  higher-order  Fast  Marching  Method  than  the  hrst- 
order  one  (see  Section  3.1.2)  is  discussed  in  Section  5.5.  Higher-order  FMM  results 
in  a  signihcant  reduction  of  errors  in  distance  estimates,  i.e.  the  diherence  between 
numerically  computed  and  true  optimal  distances.  Thus,  the  higher-order  FMM  helps 
in  dealing  with  divergence  problems  to  some  extent.  In  Section  5.6,  methods  to  deal 
with  negative  covariances  arising  due  to  the  presence  of  islands  and  due  to  the  numer¬ 
ical  error  in  computing  the  optimal  path  length  are  discussed.  These  methods,  which 
can  remove  divergence  problems,  are  discarding  the  problematic  data,  introducing  the 
process  noise  and  the  dominant  singular  value  decomposition  of  a-priori  covariance. 


5.1  Sequential  Processing  of  Observations 

A  block-diagonal  structure  of  the  observation  error  covariance  matrix  (R)  can  be 
very  advantageous  for  improving  the  computational  efficiency  and  also  for  understand- 


62 


ing  the  impact  of  the  individual  data  during  the  update  step  of  the  Kalman  hlter. 
Many  pairs  of  observations,  such  as  the  observations  from  different  instruments,  can 
have  uncorrelated  errors  resulting  in  a  block-diagonal  error  covariance  matrix  (R). 
In  fact.  Chapter  2  states  that  in  most  situations,  the  observation  matrix  for  each 
data  type  is  chosen  diagonal  with  a  uniform  non-dimensional  error  variance  e^,  i.e. 
R  =  e^I.  A  Cholesky  factorization  can  also  be  applied  to  diagonalize  the  R  matrix 
prior  to  using  sequential  processing  of  observations.  The  modihed  algorithm,  which 
takes  the  advantage  of  uncorrelated  observations  and  the  block-diagonal  structure  of 
R  is  now  discussed  (see  Parrish  and  Cohn  (1985),  Cho  et  ah  (1996)).  One  step  (one 
scale)  of  the  MSEAS  OA  equations  (Equations  2.1,  2.2,  2.3)  can  be  rewritten  as: 


K  =  Cor(x,  x)R^[iPCor(x,  x)Tr^ -|- R]  ^ 

(5.1) 

pOA  _  {I  -  KH)Cor{^,^) 

(5.2) 

=  ^  +  P[d-P^] 

(5.3) 

where  H  is  the  observation  matrix,  tjj  is  the  background  held  and  is  the  error 
covariance  of  the  estimated  held.  The  block-diagonal  matrix  R  can  be  written  as: 


R 


1 


R  = 


R2 


(5.4) 


V 


where  each  block  Rj,  j  =  1,  2, ...  J  is  apjXpj  matrix,  such  that  Pi+P2+---+Pj  =  pis  the 
total  number  of  observations.  Similarly  the  observation  vector  d  and  the  observation 
matrix  H  can  be  written  as: 


(^2 

H2 

d  = 

and  H  = 

1  j 
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Here,  dj  is  a  vector  of  size  pj  and  Hj  is  the  pj  x  n  observation  matrix  corresponding 
to  the  data  in  batch  j.  The  scalar  n  is  the  length  of  the  grid  point  position  vector  x. 
Thus,  there  are  J  batches  of  observations  such  that  the  observation  error  in  the  same 
batch  may  be  correlated  but  the  observations  in  different  batches  are  uncorrelated. 
In  order  to  sequentially  process  the  observations  by  batch,  lets  begin  by  dehning: 


Cor\x,x)Q  =  Cor-{x,x)  and  Tq  =  T.  (5.6) 

Then  the  modihed  equations  for  the  sequential  processing  algorithm  with  j  =  1,2,...N 
are: 

Kj  =  Cor{-xi,x.)j_iHj[HjCor{-K,x.)j_iHj  +  Rj]~^  (5.7) 

Cor(x,  x)j  =  {I  —  KjHj)Cor{x,x)j^i  (5.8) 

=  V’j-i +  ^i[di (5.9) 


This  gives  =  tjjj  and  =  Cor{'x,:s.)j  for  each  scale  (step)  of  the  OA  scheme. 
The  equivalence  of  Equations  5.1  -  5.3  and  Equations  5.7  -  5.9  as  proven  by  Parrish 
and  Cohn  (1985)  is  shown  in  Appendix  A.  Due  to  the  smaller  sizes  of  the  matrices 
to  invert.  Equations  5.7  -  5.9  with  the  relevant  matrix  dimensions  pj  are  signihcantly 
less  expensive  than  Equations  5.1  -  5.3  with  dimension  p. 

The  detailed  implementation  of  the  sequential  processing  Equations  (5.7  -  5.9)  is 
now  discussed.  Consider  the  case  when  all  observations  are  uncorrelated  i.e  pj  =  1, 
j=l,2,...J.  Then,  the  scalar  Rj,  row  vector  Hj  of  length  n  and  the  column  vector  Kj 
of  length  n  are  written  as: 


=  h,  =  Hj,  k,=K,  (5.10) 

The  hnal  sequential  processing  algorithm  is: 

Vj  =  Cor{x,x)j-ihj 
aj  =  hjvj  +  a^ 
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(5.11) 


V’i  =  'ipj-1  +  Pjkj 
Cor{x,x)j  =  Cor{x,x)j_i  —  kjvj 

The  computational  efficiency  can  be  further  improved  by  taking  the  sparsity  of  hj 
into  account.  For  example,  when  an  observation  falls  at  a  grid  point,  hj  will  consist  of 
a  single  one  and  the  rest  zeros.  In  this  case,  the  vj  will  be  a  column  of  Cor{x,x)j-i. 
Apart  from  increasing  the  computational  efficiency,  the  above  algorithm  also  allows 
to  estimate  the  impact  of  the  individual  observation  on  the  OA  held. 


5.2  Positive  definite  correlation  functions: 
Weiner-Khinchin  and  Bochner  theorems 

It  was  discussed  in  Chapter  2  that  the  covariance  matrix  is  generated  using  analyt¬ 
ical  correlation  functions.  Such  analytical  correlation  functions  are  termed  “positive 
dehnite  correlation  functions”  if  they  generate  positive  dehnite  covariance  matrix  us¬ 
ing  the  Euclidean  distance  for  a  simply-connected  convex  domain.  It  has  been  well 
established  using  the  Wiener-Khinchin  relationships  that  if  a  Fourier  transform  (or 
the  spectral  density  of  a  correlation  function)  is  non-negative  for  all  frequencies  then 
the  correlation  function  is  positive  dehnite  (Yaglom,  1987,  Papoulis,  1991,  Yaglom, 
2004,  Dolloh  et  ah,  2006).  The  spectral  density  (or  power  spectrum)  S{uj)  is  dehned 
as  the  Fourier  transform  of  the  correlation  function,  i.e.,  R{t)  <->■  *S'(ci;),  where: 

/OO 

R{T)e-^^^dT 

-OO 

R{t)  =  ^  r  S{oj)e^^^duj  (5.12) 

27r  J-oo 

The  above  equations  are  known  as  Wiener-Khinchin  relationships.  The  sufficient 
condition  for  the  validity  of  the  above  equation  (to  guarantee  the  existence  of  a 
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Fourier  transform  pair)  is: 


\R{T)\dT  <  oo. 


(5.13) 


The  Wiener-Khinchin  relationships  are  a  part  of  the  Wiener-Khinchin  theorem,  which 
states  that  R{t)  is  a  positive  dehnite  correlation  function  if  S{u)  >  0  V  cu.  This  is 
also  known  as  the  Bochner  theorem  (Yaglom,  1987,  2004,  Dolloff  et  ah,  2006). 


The  proof  for  the  Weiner-Khinchin  relationships  or  Bochner  theorem  (see  Papoulis 
(1991),  Yaglom  (2004),  Dolloff  et  al.  (2006))  is  described  next.  Let  T  be  the  matrix 
obtained  from  the  correlation  function  R{t).  The  proof  of  the  statement,  which  says 
that  T  will  be  a  semi-dehnite  matrix  if  the  spectral  density  5'(ci;)  is  non-negative 
>  0),  is  as  follows: 

Let  z  be  an  arbitrary  vector.  Since  (5'(ci;)  >  0)  therefore. 


z^Tz 


^  ^  tk) 

i,k 

1 

i,k 

ZTT  J —oo  ^ 


(5.14) 


Note  that  the  above  proof  assumes  that  the  correlation  function  depends  on  the 
difference  ti  —  tk  only,  which  is  the  case  for  the  Euclidean  distance  in  a  simply- 
connected  convex  domain.  But,  this  is  not  necessarily  the  case  for  the  optimal  path 
length  computed  on  multiply- connected  domains  (e.g.  using  the  Fast  Marching  or  the 
Level  Set  Method).  So,  some  computational  problems  may  be  observed  with  our  new 
OA  schemes  for  complex  multiply-connected  coastal  regions  and  archipelagos  since 
the  covariance  matrix  may  not  be  positive  dehnite.  Some  numerical  approaches  to 
deal  with  such  divergence  problems  will  be  discussed  in  this  Chapter. 
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Using  the  above  theorem,  the  necessary  condition  for  the  correlation  function 
dehned  in  Equation  2.4  to  be  positive  dehnite,  can  be  obtained.  The  conditions  are 
derived  next.  The  following  Fourier  transform  relations  (Dolloff  et  ah,  2006)  are  used: 


(5.15) 

(5.16) 


Here,  F  is  the  Fourier  transform  of  f(r).  Using  Equation  5.15,  the  following  Fourier 
transform  pairs  are  obtained: 


r.2 

.jd 

2L% 


(2„Li 

(2nLl 


1  ■ 

|2e“ 

1  ■ 
i2e' 


(5.17) 

(5.18) 


Multiplying  Equations  5.17  and  5.18  gives: 


(27rL^ 


/  2  I  2  \  r  2 


(5.19) 


Using  Equations  5.16  and  5.17,  the  Fourier  transform  of  {x^  +  y‘^)e  w.r.t.  the 
variable  x  is  given  by: 

{x^  +  y^)e  ^  ^  (27rLg)5e^~[Lg  -  +  |/^(27rLg)^e^~.  (5.20) 


Using  Equations  5.20,  5.18  and  5.16,  the  Fourier  transform  of  +  y‘^)e  w.r.t. 
the  variable  y  is  given  by: 


{x  +  y  )e 


(5.21) 


Finally,  Equations  5.21  and  5.19  are  combined  to  obtain  the  Fourier  transform  of  the 
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correlation  function  in  equation  2.4.  Thus: 


^  x^  +  y^\  - 

.  LI  )" 


{2nLi)e - ^ 


1  -  2-|  +  {ul  +  UJ2)j^ 
-^0  -^0 


(5.22) 


Thus,  Equation  5.22  shows  that  the  correlation  (Equation  2.4)  will  be  positive  dehnite 


when 


1  -  2 


Lg  +  (^^1  +^2)12 


>  0,  i.e.  when  Lq  >  \/2Le,  as  the  power  spectrum  will 
then  be  non-negative  for  all  u.  It  is  important  to  note  that  the  above  derivation 
holds  only  when  the  Euclidean  distance  is  used  along  with  the  correlation  function  in 
a  simply-connected  convex  domain. 


5.3  Divergence  issues  with  the  Kalman  filter 

Our  mapping  schemes  can  encounter  divergence  problems  associated  with  the 
Kalman  Eilter  state  and  covariance  update  equations  (Brown  and  Hwang,  1997).  If 
the  number  of  observations  is  large,  divergence  problems  can  arise  under  certain  con¬ 
ditions  due  to  truncation  errors  even  if  the  background  covariance  is  positive  dehnite. 
Like  any  numerical  procedure,  the  round  off  errors  can  lead  to  problems  as  the  num¬ 
ber  of  steps  increase.  This  can  also  occur  in  the  sequential  processing  of  observations 
(recursive  algorithm).  If  the  covariance  matrix  loses  the  property  of  positive  deh- 
niteness,  then  the  solution  of  the  Gauss  Markov  Estimation  (equivalent  here  to  the 
Kalman  Filter  update)  will  not  necessarily  be  stable.  Of  course,  such  problems  due 
to  truncation  errors  can  be  minimized  by  using  a  high-precision  arithmetic. 

Some  of  the  useful  techniques  described  by  Brown  and  Hwang  (1997)  to  counter 
the  divergence  problems  are: 

a.  Introduce  Process  Noise:  As  the  number  of  steps  increases,  the  covariance 
matrix  approaches  a  semidehnite  condition.  A  small  numerical  error  can  thus  make 
the  covariance  matrix  non-positive  dehnite,  which  may  lead  to  divergence  problems. 
A  solution  is  to  add  a  small  process  noise  (model  wise  or  in  the  OA  case  background 
noise)  to  the  diagonal  elements  of  the  covariance  matrix.  This  will  lead  to  a  degree 
of  sub-optimality  but  will  prevent  the  divergence  problems  from  arising. 


b.  Symmetrize  covariance  matrix:  The  covariance  matrix  may  become  asym¬ 
metric  due  to  roundoff  errors  and  this  asymmetry  can  grow  if  it  is  left  unchecked, 
leading  to  divergence  problems.  The  remedy  is  to  symmetrize  the  covariance  ma¬ 
trix  after  each  recursive  step.  One  approach  will  be  to  assume  the  symmetry  and  to 
use  only  the  upper  (or  lower)  triangular  part  of  the  covariance  matrix  for  all  matrix 
operations. 

c.  Use  Joseph’s  form:  If  a  large  initial  uncertainty  (covariance  matrix  with  large 
eigen  values)  is  followed  by  a  very  precise  measurement  (small  observational  error, 
R),  the  covariance  matrix  has  to  transition  from  a  large  value  to  a  small  value  in  a 
single  step  and  this  situation  can  numerically  be  very  sensitive.  The  Kalman  Filter 
update  equation  for  the  covariance  (P)  is  given  by  (see  also  Table  2.1): 

Pt  =  {I-  (5.23) 

It  is  recommended  to  use  other  covariance  update  formulas  like  the  Joseph  form 
(Grewal  and  Andrews,  1993)  which  has  a  natural  symmetry.  The  Joseph  form  is 
given  by: 


Pt  =  (I-  KtHt)Pt\t-iiI  -  KtHtf  +  KtRtKj.  (5.24) 

The  Joseph  form  (Equation  5.24)  has  a  better  numerical  behavior  (Grewal  and  An¬ 
drews,  1993)  for  unusual  numerical  situations  than  the  form  in  Equation  5.23. 
d.  U-D  factorization:  The  U-D  factorization  algorithm  (Bierman,  1977),  which  is 
mathematically  equivalent  to  the  regular  Kalman  Filter  algorithm  and  belongs  to  the 
class  of  Kalman  Filter  algorithms  known  as  square-root  hltering  (Battin,  1964,  Bier¬ 
man,  1977,  Maybeck,  1979),  has  much  better  numerical  behavior  for  a  large  number 
of  steps  (equivalent  to  the  number  of  observations).  This  algorithm,  which  will  pre¬ 
serve  the  symmetry  and  positive  dehniteness,  is  based  on  propagating  the  factors  of 
covariance  (P)  rather  than  the  full  covariance  (P).  This  algorithm  is  favored  over  the 
regular  Kalman  Filter  when  numerical  stability  is  a  concern.  The  U-D  factorization 
algorithm  is  now  described.  Note  that  it  has  some  similarity  to  the  Error  Subspace 
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Statistical  Estimation  (ESSE)  update  scheme  (see  Lermusiaux  and  Robinson  (1999), 
Lermusiaux  (1999a, b)). 


A  symmetric,  positive  definite  covariance  matrix  can  be  decomposed  into  the  form: 

P  =  UDU^.  (5.25) 

where  is  a  upper  triangular  matrix  (having  ones  along  the  major  diagonal)  and  D 
is  diagonal  (H  >  0).  The  a-priori  covariance  matrix  (Pi|o)  at  t  =  0  is  factored  using 
the  U-D  decomposition: 


-Pl|0  —  h^l|0-Dl|oh^pO- 


(5.26) 


The  Kalman  gain  (K)  and  the  state  update  are  carried  out  using  the  usual  Kalman 
Filter  equations  by  replacing  the  a-priori  covariance  using  Equation  5.26.  The  covari¬ 
ance  update  form  is: 


Pt  =  Pt\t-i-Pt\t-iH{ty{H{t)Pt\t-iH{ty  +R{t))-^H{t)Pt\t-i 


—  Pt\t-1 


Ptit-iH{t)^H{t)Ptit-, 


a 


(5.27) 


where  the  scalar  a  is  given  by: 


a  =  H{t)Pt\t-iH{tf  +  R{t) 


(5.28) 


Note  that  in  the  sequential  processing  of  observations,  a  single  observation  is  taken  in 
each  step  of  processing  since  the  observational  error  matrix  (i?)  is  diagonal.  Therefore, 
a  is  a  scalar.  In  the  general  case,  when  R  is  a  block  diagonal  matrix,  the  covariance 
update  form  is: 


Pt  —  Pt\t-i  —  Pt\t-iH{fy"(^  ^p[{t)Pt\t-i  (5.29) 
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The  covariances  P*  and  Pt\t-i  can  be  replaced  by  their  U-D  factors.  This  gives: 


UtD,U,^  = 


a 


a 


T\T-[ 


t\t-l 


(5.30) 


The  bracketed  term  in  Equation  5.30  is  symmetric,  so  it  can  be  factored  using  U-D 
factorization.  Thus: 


a 


U,D,U,^. 


(5.31) 


This  gives: 


UtDtU^  = 

=  {Ut\t-iUt)Dt{Ut\t.,Utf  (5.32) 

Since  {Ut\t-iUt)  is  an  upper  triangular  and  Dt  is  diagonal,  therefore: 

Ut  =  Utit-iUt 

Dt  =  Df  (5.33) 

These  equations  are  the  same  as  the  ESSE  update  equations  (see  Lermusiaux  and 
Robinson  (1999),  Lermusiaux  (1999a,b)). 


For  Objective  Analysis,  using  the  sequential  processing  of  observations,  the  co- 
variance  update  in  terms  of  U-D  factors  Uj  and  Dj  (j=l,2,...J,  J  is  the  number  of 
observations)  is: 


Dj+i  =  Dj. 
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(5.34) 

(5.35) 


where,  Uq  and  Dq  can  be  expressed  in  terms  of  the  background  covariance  (Cor(x,x)): 


UqDqUq  =Cor(x,  x)  (5.36) 

As  mentioned,  the  U-D  factorization  helps  in  dealing  with  divergence  issues  by  ensur¬ 
ing  that  the  covariance  matrix  always  remains  positive  dehnite.  In  situations  where 
there  are  adequate  process  noise  entering  the  system  states,  the  usual  Kalman  Filter 
update  equation  (Equation  5.23)  is  also  capable  of  handling  all  the  divergence  issues 
and  the  algorithm  based  on  U-D  factorization  may  not  be  necessary. 


5.4  Comparison  of  Computational  Costs 

For  a  2-dimensional  domain  with  N  points  in  each  direction,  a  comparison  of  the 
operation  count  for  computing  the  optimal  distance  from  a  data  location  to  all  other 
grid  points  in  the  domain  using  different  Methods  is  given  in  Table  5.1. 


Method 

Operation  Count 

Level  Set  Method 

0(iV3) 

Fast  Marching  Method 

0{NHogN) 

Dijkstra’s  Method 

0(iV3) 

Table  5.1:  Comparison  of  the  operation  count  for  the  optimal  distance  obtained  using 
LSM,  FMM  and  Dijkstra’s  Method. 


There  are  a  total  of  grid  points  at  each  level  and  the  operation  count  for  LSM 
is  obtained  from  an  optimistic  guess  that  LSM  will  take  roughly  N  steps  to  converge. 
In  reality,  the  iterations  can  take  much  longer  to  converge,  and  therefore  LSM  is  not 
a  very  efficient  method  for  computing  the  optimal  distance  to  perform  OA.  On  the 
other  hand,  FMM  is  an  efficient  technique  which  requires  a  fast  method  to  locate 
the  smallest  value  grid  point  in  the  narrow  band.  The  Min-Heap  data  structure  with 
backpointers  (Sedgewick,  1988)  is  employed  to  efficiently  locate  the  grid  point  with 
the  minimum  value.  The  total  work  done  in  the  DownHeap  and  UpHeap  operations, 
which  ensure  that  the  updated  quantities  do  not  violate  the  heap  properties,  is  0(log 
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N).  Thus  2-dimensional  FMM  with  N  grid  points  in  each  direction  has  an  operation 
count  of  N'^logN,  which  is  a  signihcant  improvement  over  LSM.  It  has  also  been 
observed  that  FMM  requires  less  computational  time  (approximately  15  %)  than  the 
SDE  approach  proposed  by  Lynch  and  McGillicuddy  (2001).  Thus,  the  OA  based 
on  FMM  is  computationally  the  most  efficient  technique  for  mapping  in  multiply- 
connected  domains. 


5.5  Higher  order  Fast  Marching  Method 


In  a  domain  with  no  islands  or  landforms,  the  optimal  path  length  obtained  using 
the  FMM/LSM  should  ideally  be  equal  to  the  Euclidean  distance.  But  the  numer¬ 
ical  estimation  of  the  optimal  path  length  using  the  EMM/LSM  has  discretization 
errors  and  this  leads  to  an  inaccurate  estimation  of  the  optimal  path  length.  We  have 
discussed  that  the  Weiner  Khinchin  and  Bochner  theorems  are  valid  for  covariances 
computed  using  the  Euclidean  distance  in  a  simply-connected  convex  domain.  The 
covariance  matrix  may  no  longer  be  positive  dehnite  due  to  the  inaccurate  computa¬ 
tion  of  the  optimal  path  length  by  FMM/LSM  or  due  to  the  presence  of  islands.  This 
may  lead  to  divergence  problems  in  the  resultant  field  maps  based  on  the  FMM/LSM 
scheme  as  shown  in  Figure  C-38.  Specihcally,  the  question  which  motivates  this  Sec¬ 
tion  is:  What  are  the  computational  errors  in  optimal  path  lengths  computed  using 
the  FMM/LSM  and  how  can  they  be  reduced?  In  this  section,  we  introduce  the 
higher  order  Fast  Marching  Method  which  will  reduce  errors  in  the  estimation  of  the 
optimal  path  length. 

The  Fast  Marching  Method  presented  in  Section  3.2  is  a  hrst  order  scheme,  since 
the  hrst  order  discretization  form  (Equation  3.10)  of  the  Eikonal  equation  (Equation 
3.7)  was  used.  A  different  implementation  of  FMM  with  higher  accuracy  (Sethian, 
1999a,b)  is  discussed  here.  Note  that  the  second  order  backward  approximation  to 
the  hrst  derivative  is  given  by: 


T 

±  rp 


3Ti  —  4Tj_i  -|-  Tj_2 
2Ax 


D~ 


ICrjn  _|_  ^  rj  —  X  —  Xrp 

2 


(5.37) 
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Similarly,  the  second  order  forward  approximation  to  the  first  derivative  is  given 
by: 


T 

■J-  T. 


‘iTi  —  4Tj+i  +  Ti 


U+2 


2Ax 


Ax 


JJ+X+Xrjn 


(5.38) 


Here  D  ^  and  are  the  hrst  order  forward  and  backward  approximations  for  the 
hrst  derivative,  respectively  (Equation  3.6),  D~^~^  =  _D+®+^  =  . 

Consider  the  switch  functions  dehned  by: 


switch,  A 


switch^'' 


1  if  Tj_2j  and  Tj_ij  are  known  (‘Alive’)  and  Tj_  2j  <  Ti—lj 
0  otherwise 

1  if  Tj+2j  and  Tj+ij  are  known  (‘Alive’)  and  Ti+2j-  <  Tj+ij 
0  otherwise 


(5.39) 


Similar  functions  are  dehned  in  the  y  direction.  The  higher  accuracy  scheme  attempts 
to  use  a  second  order  approximation  for  the  derivative  whenever  the  points  are  tagged 
as  ‘alive’  (the  points  inside  the  band  where  the  value  of  the  arrival  time  function  is 
frozen:  see  Section  3.2)  but  reverts  to  the  hrst  order  scheme  otherwise. 

The  modihed  discretization  equation  for  the  higher  accuracy  EMM  is  given  by: 


^  max{[D;^^T  +  switch^^^D^^-^T],-[D±^T  -  switch±^ ^D±^+^T],0f  ^ 

+ 

^  max{[D^yT  +  switch;^ ^D-/-yT], -[D+yT  -  switch^^ ^D^y^yT],0)^  ^ 


1 


F.2 


(5.40) 


It  should  be  noted  that  the  above  scheme  is  not  necessarily  a  second  order  scheme. 
The  accuracy  of  the  above  scheme  depends  on  how  often  the  switches  evaluate  to 
zero  and  how  the  number  of  points  where  the  hrst  order  method  is  applied  changes 
as  the  mesh  is  rehned.  When  the  number  of  points  where  the  hrst  order  method 
is  applied  is  relatively  small  (occurs  only  near  the  coastlines),  the  error  is  reduced 
considerably  by  using  the  higher  accuracy  EMM.  This  is  clearly  evident  from  Figures 
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C-39  and  C-40,  which  show  a  comparison  of  the  distances  and  correlation  (difference 
and  normalized  difference  from  Euclidean  distance)  obtained  using  the  hrst  order  and 
higher  accuracy  EMM  schemes  in  the  trivial  case  when  a  scalar  speed  function  F  is 
equal  to  one  everywhere  in  the  domain.  It  should  also  be  noted  that  a  third  or  higher- 
order  approximations  for  the  derivative  can  similarly  be  used  to  construct  more 
accurate  EMM  schemes,  but  this  will  increase  the  computational  cost.  Figures  C-39 
and  C-40  also  show  that  the  relative  error  in  the  optimal  distance  computed  using 
EMM  is  higher  near  the  data  point  and  it  reduces  as  the  distance  increases.  To  keep 
the  computational  cost  low  and  a  uniform  relative  error,  one  can  use  higher  accuracy 
EMM  near  the  data  point  and  then  progressively  shift  to  the  lower  order  schemes  as 
the  distance  increases. 

The  higher  order  Fast  Marching  Method  has  been  used  to  minimize  errors  in 
the  estimation  of  the  optimal  path  length  in  Philippines  Archipelago.  Figure  C-41 
clearly  shows  that  the  use  of  higher  order  Fast  Marching  Method  has  attenuated  the 
divergence  problems  compared  to  the  hrst  order  EMM.  The  divergence  problems  do 
not  vanish  completely  because  of  the  presence  of  landforms  and  due  to  discretization 
errors  associated  with  higher  order  EMM.  We  introduce  other  methods  to  deal  with 
such  divergence  problems  for  multiply-connected  coastal  domains  in  Section  5.6. 


5.6  Positive  Definite  covariance  matrix  for  com¬ 
plex  multiply-connected  coastal  regions 

We  discussed  in  Section  5.2  that  the  Weiner-Khinchin  and  Bochner  theorems 
are  valid  for  the  background  error  covariance  matrix  computed  using  the  Euclidean 
distance  for  simply-connected  convex  domains.  The  matrix  may  become  negative  due 
to  the  discretization  error  in  the  optimal  path  length  computed  using  FMM/LSM.  The 
use  of  higher  order  EMM  was  therefore  recommended  in  Section  5.5.  The  covariance 
matrix  may  also  become  negative  due  to  the  presence  of  islands  and  coastlines. 

Specihcally,  the  question  which  motivates  this  Section  is:  What  are  the  compu- 
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tational  issues  including  non-positive  definite  covariances  that  arise  in  a  multiply- 
connected  coastal  domain  and  how  can  they  be  remedied?  The  presence  of  islands 
and  archipelagos  results  in  stretching  of  the  Euclidean  path  in  the  physical  domain 
which  can  potentially  make  the  covariance  matrix  negative.  A  negative  covariance 
matrix  results  in  divergence  problems  associated  with  the  Kalman  Filter  (which  is 
linked  to  the  algorithm  for  sequential  processing  of  observations).  Examples  of  this 
will  be  shown  next.  Possible  remedies  are  then  discussed  starting  with  curvilinear 
coordinates  (Section  5.6.1).  Then,  other  methods  to  deal  with  the  negative  covari¬ 
ance  matrix  due  to  the  presence  of  islands  and  due  to  the  inaccurate  estimation  of 
the  optimal  path  length  using  the  EMM  or  LSM  are  discussed  (Section  5.6.2). 

Consider  the  idealized  multiply-connected  domain  having  an  island,  shown  in 
Figure  C-42.  This  domain  has  12  grid  points  which  are  marked  as  ocean  points 
and  4  grid  points  which  are  marked  as  land  points.  The  length  of  the  optimal  path  is 
computed  exactly  to  form  the  covariance  matrix  to  keep  it  untouched  from  effects  due 
to  discretization  errors  in  the  FMM/LSM.  The  positive-dehnite  correlation  function 
Cor(r)  =  exp  with  L=2  is  used  to  form  the  covariance  matrix.  We  hnd  that 

the  covariance  matrix  is  not  positive  dehnite.  The  maximum  eigen  value  for  the 
covariance  matrix  is  6.3345  while  the  minimum  is  -0.0504.  This  idealized  example 
clearly  shows  that  the  covariance  matrix  based  on  the  optimal  path  length  for  a 
complex  multiply-connected  region  may  not  necessarily  be  positive  dehnite. 

5.6.1  Curvilinear  grids 

If  a  very  high  accuracy  EMM  is  used,  the  optimal  distance  computed  using  this 
high  accuracy  EMM  is  equivalent  to  the  Euclidean  distance  only  for  simply-connected 
convex  domain.  Thus,  to  ensure  that  the  Weiner-Khinchin  relationship  or  Bochner 
theorem  hold,  a  mapping  technique  based  on  using  curvilinear  coordinates  can  be 
used  (Cebeci  et  ah,  2005).  However,  this  will  stretch  the  correlation  scales  as  we 
will  describe  now.  First,  an  irregular  region  having  islands  in  the  physical  plane  in 
the  Cartesian  (x,y)  or  the  polar  (r,0)  coordinates  is  mapped  on  the  curvilinear  {(,  p) 
coordinates  such  that  the  mapped  region  is  a  simply-connected  convex  domain. 
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A  simple  example  of  such  a  mapping  in  a  domain  having  a  circular  island  (which 
is  mapped  to  a  rectangular  domain  and  islands  are  absent  in  the  mapped  region)  is 
illustrated  in  Figure  C-43.  The  domain  having  a  circular  island  in  the  polar  (r,0) 
coordinates  is  mapped  on  the  curvilinear  coordinates.  The  optimal  distance 

in  this  transformed  coordinate  system  will  be  the  Euclidean  distance  for  which  the 
Weiner-Khinchin  and  Bochner  theorems  hold.  These  optimal  distances  also  satisfy 
coastline  constraints  i.e.  there  is  no  direct  relationship  across  landforms.  For  example, 
the  distance  between  the  locations  ((0,7]3)  and  (CO,?]?)  (which  are  separated  by  the 
island  in  the  polar  coordinates)  in  the  curvilinear  coordinates  (Figure  C-43)  will  be 
equal  to  the  length  of  the  optimal  path  in  the  polar  coordinates  which  is  different 
from  the  across-land  Euclidean  distance  between  the  locations  (rO,03)  and  (rO,07)  in 
the  polar  coordinates.  It  should  also  be  noted  that  the  distance  between  the  locations 
(C0,?73)  and  ((0,?]7)  is  equal  to  the  distance  between  the  location  (C3,?]3)  and  (C3,?]7) 
in  the  curvilinear  coordinates.  This  is  not  the  case  in  the  polar  coordinates  where 
the  length  of  the  optimal  path  between  (rO,03)  and  (rO,07)  is  smaller  than  the  length 
of  the  optimal  path  between  (r3,6'3)  and  (r3,6*7).  This  corresponds  to  decreasing  the 
grid  resolution  away  from  the  island.  Thus,  the  mapping  in  transformed  coordinates 
will  lead  to  deformations  of  the  length  scales  specihed  for  the  Objective  Analysis  in 
the  physical  region.  Also,  the  mapping  becomes  complicated  when  the  domain  has 
multiple  islands  of  distorted  shapes.  Therefore,  alternative  methods  are  reqnired  to 
remove  problems  arising  dne  to  the  presence  of  islands. 

5.6.2  Other  methods 

Other  methods  which  are  very  nsefnl  in  removing  the  divergence  problems  (Figure 
C-44  (Top-Left))  due  to  the  negative  covariance  matrix  are: 

a.  Discarding  the  problematic  data:  One  method  to  deal  with  the  problem  of  a 
negative  covariance  matrix  is  to  discard  the  data  corresponding  to  the  negative  valnes 
of  i7jC'or(x,  x)j_ii7j.  Even  thongh,  this  will  ensnre  that  there  are  no  divergence 
issnes  in  the  resnltant  OA,  this  method  is  a  not  the  most  acceptable  one  since  the 
information  in  the  data  is  discarded  entirely.  The  held  map  obtained  by  discarding 
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the  problematic  data  is  shown  in  Figure  C-44  (Top- Right).  Clearly,  the  divergence 
problems  are  removed  but  loosing  all  the  information  in  the  data  is  not  acceptable. 

b.  Introducing  process  noise:  As  discussed  in  Section  5.3,  adding  a  small  process 
noise  to  the  diagonal  elements  of  the  covariance  matrix  helps  in  dealing  with  the 
divergence  problems  associated  with  a  negative  covariance  matrix.  The  disadvantage 
is  that  the  process  noise  introduced  will  lead  to  a  degree  of  sub-optimality.  It  is  often 
a  more  acceptable  method  compared  to  discarding  the  data.  Once  again,  the  held 
map  obtained  by  introducing  the  process  noise  is  free  from  divergence  problems  and 
the  plot  is  shown  in  Figure  C-44  (Bottom-Left). 

c.  Dominant  Singular  Value  Decomposition  (SVD)  of  a-priori  covariance: 

To  construct  the  held  maps  using  the  OA  techniques,  the  full  covariance  matrix  is 
not  required.  The  computation  of  the  full  covariance  matrix  (Cor(x,x))  is  expensive 
and  it  is  therefore  rarely  done  for  the  OA  in  complex  coastal  regions.  The  necessary 
requirement  to  obtain  the  held  maps  is  the  matrix  corresponding  to  the  grid  and  the 
data  point  covariance  (Cor(x,X)).  The  divergence  problems  in  the  Kalman  update 
or  in  the  sequential  processing  of  observations  can  be  removed  by  hrst  obtaining 
the  singular  value  decomposition  (SVD)  of  Cor(x,X)  and  then  reconstructing  the 
new  covariance  matrix  by  retaining  only  the  dominant  singular  values  and  setting 
the  smaller  singular  values  (less  than  1  percent  of  the  maximum  singular  value)  to 
zero.  The  above  procedure  will  make  the  covariance  positive  dehnite.  It  has  been 
verihed  that  the  magnitude  of  the  negative  eigen  values  in  the  covariance  matrix  is 
very  small  compared  to  the  magnitude  of  the  maximum  eigen  value.  This  verihcation 
establishes  that  the  use  of  the  dominant  singular  value  decomposition  method  is  the 
most  acceptable  method  to  remove  the  divergence  problems  in  the  update  because  it 
looses  the  least  information  contained  in  the  data.  Once  again,  the  held  map  obtained 
by  dominant  singular  value  decomposition  (SVD)  of  a-priori  covariance  is  free  from 
divergence  problems  and  the  plot  is  shown  in  Figure  C-44  (Bottom- Right).  It  has  also 
been  verihed  that  helds  obtained  by  introducing  the  process  noise  and  helds  obtained 
by  dominant  singular  value  decomposition  of  the  a-priori  covariance  are  very  similar. 

Another  important  and  challenging  problem  in  oceanography  is  estimating  the 
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scales  from  the  data.  New  methodologies  for  scale  estimation  based  on  the  structure 
function  method,  short  term  Fourier  transform  (STFT)  and  second  generation  wavelet 
(SGW)  analysis  are  discussed  in  the  next  Chapter. 
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Chapter  6 


Adaptive  Scale  Estimation 


A  major  but  apparently  simple  challenge  in  geophysical  fluid  dynamics  at  sea 
is  the  estimation  of  the  dominant  spatial-time  scales  in  held  measurements.  This 
may  seem  basic  for  scientists  and  engineers  not  used  to  ocean  data,  but  due  to  the 
multi-scale,  turbulent  and/or  intermittent  ocean  dynamics  and  due  to  the  irregular 
properties  of  the  data,  it  is  a  very  common  question  in  ocean  meetings.  In  this 
chapter,  we  investigate  new  methods  for  adaptive  spatial-time  scale  estimation  from 
irregular  ocean  held  measurements.  Solving  this  question  would  significantly  help  in 
better  understanding  and  sampling  of  ocean  processes. 

New  adaptive  schemes  to  learn  the  largest  and  the  most  energetic  scales  based 
on  structure  functions  (Denman  and  Freeland,  1985)  and  on  non-linear  least  square 
fitting  methods  (Appendix  B)  are  first  outlined  in  Section  6.1.  An  example  of  an 
adaptive  scheme  based  on  short  term  Fourier  transforms  (STFT),  which  can  be  used 
for  non-stationary  signals,  is  illustrated  in  Section  6.2.  It  is  applied  to  an  idealized 
dataset  generated  using  a  chirp  signal,  in  which  the  frequency  and  the  wave  number 
varies  with  time.  In  Section  6.3,  we  illustrate  another  new  adaptive  scale-estimation 
scheme  based  on  second  generation  wavelets  (Sweldens,  1998,  Jansen  and  Oonincx, 
2005).  Second  generation  wavelets  are  applicable  to  both  irregularly  sampled  and 
non-dyadic  data  sets.  This  is  important  since  we  would  like  to  estimate  scales  in 
the  data  without  mapping  the  data.  To  map  the  data,  one  needs  scale  estimates,  as 
described  in  previous  chapters.  Ultimately,  the  goal  of  all  new  methods  investigated 
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below  would  be  to  create  a  map  of  time  and  space  scales  that  evolves  as  ocean  data 
are  collected  or  are  fed  to  the  scale  estimation  scheme,  all  without  requiring  to  map 
the  data. 

6.1  A  new  adaptive  scale  estimation  method  using 
structure  functions 

The  new  adaptive  algorithm  derived  in  this  section  can  be  employed  to  estimate 
the  spatial-time  scales  from  the  available  oceanic  data  held.  The  approach  proposed 
by  Denman  and  Freeland  (1985)  is  utilized  to  obtain  the  isotropic  structure  function 
{B{r))  from  the  data  held  {6).  The  structure  function  {B{r))  is  dehned  by: 

B,{r)  =  E[{e{xi  +  r)  -  e{xi)f]  +  2N 

=  E[9{xi  +  r)9{xi  +  r)]  +  E[9{xi)9{xi)]  -  2E[9{xi  +  r)9{xi)\  +  2N 
=  2[F(0) -F(r) +  iV] 

=  2iV +  2F(0)[1 -C'or(r)]  (6.1) 

where  E{r)  =  E[9{xi  +  r)9[xi)],  Cor{r)  =  E{r)/E{0)  and  N  is  the  noise  variance. 
We  have  assumed  that  errors  in  measurements  9{xi  +  r)  and  9{xi)  are  independent. 
So  the  variance  of  the  error  in  {9{xi  +  r)  —  9{xi))  is  given  by  2N.  The  noise-free 
correlation  fnnction  with  properties  Cor{0)  =  1  and  Cor{r)  — 0  as  r  — cx)  has  been 
given  in  Equation  2.4. 

Onr  new  adaptive  learning  algorithm  utilizes  the  structure  function  B  (r)  obtained 
from  the  available  data  to  learn  the  nnknown  scale  parameters  by  using  the  non¬ 
linear  least  sqnares  htting  method  (Appendix  B).  An  adaptive  learning  algorithm  for 
estimating  the  length  scale  (L)  and  the  time  scale  (r)  is  as  follows: 

1.  Assnme  an  initial  approximation  for  the  scales.  Consider  the  data  in  the  time 
window  At  =  Ts/A 

2.  Obtain  the  structure  function  from  the  data  available  in  the  time  interval  At 
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and  the  spatial  window  \x  —  Xi\  <  for  all  the  grid  points. 

3.  Length  scales  (Lf  and  Lf)  are  obtained  by  the  non-linear  least  square  £t  of  the 
structure  function  from  the  data  to  the  analytical  structure  function. 

4.  The  following  learning  equation  is  used  to  obtain  the  new  length  scales. 

=  Lf  +  ;r(Lf  -  Lf ) 

Li‘+‘'  =  Lf  +  /r(Lf  -  Lf)  (6.2) 

Here,  Ir  is  the  learning  rate. 

5.  A  similar  analysis  is  repeated  for  obtaining  the  time  scales. 

It  should  be  noted  that  more  than  one  spatial  or  time  scales  can  be  estimated  by  ap¬ 
propriately  selecting  the  spatial  and  temporal  window.  The  data  in  the  corresponding 
spatial-temporal  window  can  be  used  for  estimating  the  large  and  small  scales. 

Scales  estimated  from  the  adaptive  learning  algorithm  based  on  the  structure  func¬ 
tion  in  the  complex  domain  of  Philippines  Archipelago  using  the  Melville  (Summer 
2007)  exploratory  cruise  data  (Figure  C-23)  are  shown  in  Figure  C-45.  Small  scale 
estimates  on  26*^  June,  2007  are  obtained  from  the  exploratory  cruise  data  collected 
on  and  before  26*^  June,  2007.  Scale  estimates  are  available  in  the  Bohol  Sea  and 
near  the  Panay  island.  As  more  data  gets  collected,  the  new  scale  estimate  plot  on 
15*^  July,  2007  shows  that  estimates  are  also  available  in  the  Mindoro  strait.  It  can 
also  be  observed  that  scales  in  the  Bohol  Sea  are  slightly  smaller  due  to  the  change 
in  wind  circulation  during  the  period  between  26*^  June,  2007  and  15*^  July,  2007. 

6.2  A  new  adaptive  scale  estimation  method  using 
short  term  Fourier  transforms  (STFT) 

Another  adaptive  learning  scheme  based  on  short  term  Fourier  transforms  (STFT) 
is  proposed  in  this  section  to  estimate  the  spatial-time  scales,  but  assuming  a  regularly 
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spaced  and  idealized  dataset.  This  allows  us  to  evaluate  capabilities  of  the  STFT 
method. 

Idealized  data  corresponding  to  the  chirp  signals  (swept  sinusoids)  has  been  gen¬ 
erated  and  the  adaptive  algorithm  is  employed  to  learn  the  instantaneous  frequency 
of  the  chirp  signal.  The  data  has  been  generated  at  N  [N  =  40000)  time  instances 
for  the  spatial  domain  (0  <  a;  <  X  =  50)  and  the  time  domain  (0  <  t  <  T  =  50). 
The  data  function  is  given  by: 


u{xi,  ti)  =  ?,sin{ki{xi)xi  +  Ui{ti)ti}  +  sin{k-2{xi)xi  +  uj2{ti)ti} 


where, 


ti  =  T 


i  -  1 

~w 


ki{xi) 


271 


Xi 


271  271  . 


1.5  ^2X^.5  1.5^ 

27r  ti  .271  27r , 

2.5  2TX.5  2.5^ 


Xi  =  X 


mod{i  —  1,  \/N) 
^/N-l 


271  Xi  271  271 
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2'.,/  6  12' 


(6.3) 


The  instantaneous  wave  numbers  (kunst,  k2inst),  frequencies  (uJunst,  ^2inst),  length 
scales  {Liinst,  L2inst)  and  time  scales  {Tunsti  T2inst)  corresponding  to  this  chirp  signal 
are  given  by: 
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klinstyXi)  =  T  )^^2~5  T~5  ’  k2inst{Xi) 
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The  stepwise  algorithm  based  on  STFT,  which  has  been  implemented  to  adaptively 
learn  the  length  and  time  scales,  is  as  follows: 


1.  Assume  an  initial  approximation  for  the  small  and  large  spatial  scales  {Lg,  Li) 
and  the  small  and  large  time  scales  (T*,  T/)  at  all  the  grid  points. 
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2.  Consider  the  data  in  a  time  interval  which  is  significantly  smaller  than  the  small¬ 
est  time  scale  {Tg/S  >  At).  Apply  the  short  term  Fourier  transform  (STFT)  in 
the  spatial  window  {\x  —  Xi\  <  Li)  and  analyze  the  corresponding  wave  number 
spectrum.  The  two  wave  numbers  corresponding  to  the  peak  amplitudes  of  the 
wave  number  spectrum  will  give  an  estimate  of  the  small  (Tf )  and  the  large 
(Lf )  spatial  scales  in  the  signal.  The  new  length  scales  are  obtained  using  the 
learning  equations: 


Lf+'l  =  if  +  lr(Lf  -  Lf) 

Li‘+‘'  =  11  +  lr(Ll  -  Lj)  (6.5) 

Here,  Ir  is  the  learning  rate. 

3.  A  similar  analysis  is  repeated  for  obtaining  the  time  scales.  STFT  is  applied 
in  the  time  window  {\t  —  ti\  <  Ti)  to  analyze  the  data  frequency.  The  two 
frequencies  corresponding  to  the  peak  amplitudes  of  the  frequency  spectrum 
will  give  an  estimate  of  the  small  {Tf)  and  the  large  {Tf)  time  scales  in  the 
signal.  The  new  time  scales  can  be  obtained  by  using  the  learning  equations 
(similar  to  Equation  6.5).  The  process  is  repeated  to  adaptively  learn  the  new 
length  and  time  scales  as  the  new  data  is  received  for  the  next  time  interval  At 
{Tg/8  >  At). 

Plots  for  the  small  and  the  large  length  scales  for  chirp  signal  data  (Equation  6.3) 
computed  using  the  adaptive  learning  algorithm  described  above,  are  shown  in  Figure 
C-46.  A  comparison  has  been  made  between  the  analytical  instantaneous  length  scales 
(Equation  6.4)  and  the  scales  adaptively  learned  from  the  data  with  a  learning  rate, 
Ir  =  0.1.  The  accuracy  of  scales  depends  on  the  resolution  of  Fourier  spectrum. 
The  analytical  small  scales  and  the  small  scales  obtained  using  the  adaptive  learning 
algorithm  compare  well,  since  the  STFT  wave  number  resolution  is  good  (see  Figure 
C-46  (Top-Left)).  Jumps  observed  in  the  small  scales  are  due  to  the  finite  resolution 
of  wave  numbers.  However,  for  the  large  scales,  the  STFT  wave  number  resolution  is 
not  as  good  as  the  small  scales,  and  the  scale  estimate  obtained  using  the  adaptive 
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learning  algorithm  does  not  match  well  with  the  analytical  instantaneous  scales  (see 
Figure  C-46  (Top-Right)).  Similar  trends  are  observed  from  plots  for  time  scales 
shown  in  Figure  C-46.  There  is  a  good  comparison  between  plots  for  small  time-scales 
obtained  using  the  learning  algorithm  with  a  learning  rate  0.1  and  the  instantaneous 
small  time  scales  corresponding  to  the  chirp  signal  data  (see  Figure  C-46  (Bottom- 
Left)).  But,  the  large  time  scales  do  not  compare  well  since  the  frequency  resolution 
is  not  good  (see  Figure  C-46  (Bottom-Right)). 

Results  obtained  above  are  encouraging  and  a  possible  approach  to  improve  the 
scale  estimates  is  by  adaptively  learning  scales  using  the  wavelet  analysis  (since  the 
signal  is  not  stationary  and  the  frequency  changes  with  time).  To  understand  effects 
of  varying  the  learning  rate  (Ir),  spatial-time  scales  have  been  obtained  for  the  same 
chirp  signal  data  with  a  learning  rate  Ir  =  0.2.  With  an  enhanced  learning  rate,  the 
scales  are  now  more  sensitive  to  the  incoming  data.  The  plots  for  the  large  spatial 
scales  and  the  large  time  scales  (see  Figure  C-47  (Right  panel))  obtained  using  the 
adaptive  learning  algorithm  with  learning  rate  Ir  =  0.2  show  that  the  oscillations  in 
the  estimated  length  and  time  scales  are  severe  compared  to  the  adaptive  learning 
with  learning  rate  Ir  =  0.1.  The  advantage  of  having  a  larger  learning  rate  is  that 
scales  will  converge  to  the  true  scales  in  a  small  number  of  steps,  but  at  the  cost  of 
higher  sensitivity  to  the  incoming  data. 


Another  data  set  has  been  generated  at  N  (iV  =  40000)  time  instances  for  the 
spatial  domain  (0  <  a;  <  X  =  50)  and  the  time  domain  (0  <  t  <  T  =  50)  to 
validate  the  adaptive  algorithm  based  on  STFT  for  learning  the  spatial-time  scales. 
The  wavenumber  and  frequency  for  the  data  function  (Data  Set  2)  are  given  by: 
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The  true  spatial-time  scales  are  given  by: 
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Our  new  STFT-based  adaptive  algorithm  for  learning  the  spatial-time  scales  is  em¬ 
ployed,  and  the  estimated  scales  are  compared  with  the  true  scales.  Plots  comparing 
the  spatial-time  scales  obtained  using  the  learning  algorithm  and  the  true  spatial-time 
scales  are  shown  in  Figure  C-48.  The  estimated  small  spatial-temporal  scales  com¬ 
pare  well  with  the  true  small  spatial-temporal  scales  (see  Figure  C-48  (Left  panel)) 
but  oscillations  are  observed  in  the  large  spatial-temporal  scale  estimates  (see  Figure 
C-48  (Right  panel))  due  to  limited  resolution  by  STFT.  These  results  indicate  that 
the  learning  algorithm  suggested  in  this  section  is  a  reasonable  approach  to  estimate 
the  spatial-temporal  scales.  Further  improvement  in  accuracy  of  the  estimated  scales 
can  be  expected  by  using  wavelet  analysis. 

Note  that  STFT  can  only  be  used  on  data  signals  which  are  regularly  spaced.  But 
the  oceanic  data  is  collected  by  platforms  whose  positions  vary  irregularly  with  time. 
Therefore,  the  learning  algorithm  based  on  the  structure  function  (Section  6.1)  and 
the  algorithm  based  on  second  generation  wavelets  (Section  6.3)  will  be  more  relevant 
to  the  ocean  data.  The  above  work  was  mainly  carried  out  to  quantify  capabilities  of 
the  STFT  method. 


6.3  A  new  adaptive  scale  estimation  method  using 
second  generation  wavelets 

Scales  in  the  ocean  vary  with  time.  The  scale  estimation  is  thus  a  non-stationary 
problem,  which  makes  the  use  of  Fourier  Transform  less  adequate.  Short  term  Fourier 
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transforms  (STFT)  and  first  generation  wavelets  are  applicable  to  non-st  at  ionary 
signals  and  are  proven  to  be  very  useful  in  engineering  and  computer  science,  but 
they  require  regularly  spaced  dyadic  data  which  is  not  the  case  with  ocean  data. 
The  ocean  data  is  irregular  and  non-dyadic.  The  use  of  a  regular  held  obtained  by 
mapping  the  data  should  also  be  avoided  if  possible,  since  the  mapping  procedure 
introduces  artihcial  scales  and  therefore  scales  estimated  from  mapped  helds  will  be 
different  from  scales  estimated  using  the  ocean  data. 

Wavelets  are  mathematical  functions  that  divide  the  data  into  different  frequency 
components  and  each  component  is  then  studied  with  a  resolution  that  matches  its 
scale.  Two  well-known  admissible  wavelets  (Figure  C-49)  which  have  been  used  in 
several  applications  are: 

a.  Haar  Wavelet:  The  Haar  wavelet  is  dehned  by: 


%Ij{x)  =  { 
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a.  Mexican  Hat  Wavelet:  Mexican  hat  wavelet  is  the  second  derivative  of  the 
Gaussian  function.  It  is  given  by: 
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Second  generation  wavelet  (SGW)  decomposition  (Sweldens,  1998,  Jansen  and  Oon- 
incx,  2005),  which  is  evaluated  directly  on  irregular  locations  and  does  not  need  dyadic 
data,  is  now  discussed.  The  approach  can  also  automatically  adapt  to  the  data  on 
hnite  intervals,  without  requiring  artihcial  techniques  like  extension  of  data  near  the 
boundaries.  SGW  is  based  on  using  the  so  called  lifting  scheme.  Lifting  can  be  viewed 
as  a  method  to  construct  a  wavelet  transform,  which  has  the  ability  to  enhance  the 
hrst  generation  wavelet  transform,  by  adding  desirable  properties.  A  description  of 
the  lifted  version  of  the  Haar  wavelet  transform  (Sweldens,  1998,  Jansen  and  Oonincx, 
2005)  is  now  provided. 


A  Haar  decomposition  takes  the  fine  scale  scaling  coefficients  as  inputs.  It  can  be 
assumed  that  these  fine  scale  coefficients  are  observations.  The  resolution  level  cor¬ 
responding  to  observations  is  the  finest.  Second  generation  wavelet  analysis  proceeds 
in  successive  steps  by  extracting  the  coarser  scale  trends  in  each  successive  step.  A 
lifting  scheme  consists  of  three  operations:  a  split,  followed  by  a  sequence  of  dual 
and  primal  lifting  operations.  As  shown  in  Figure  C-50,  the  splitting  step  divides 
the  input  signal  /  into  two  disjoint  sets:  observations  with  odd  indices  (71)  and  even 
indices(Ai).  Each  step  of  the  Haar  decomposition  computes  averages  and  differences 
of  the  adjacent  input  values.  Let  Sj+i^k  be  the  input  at  scale  j  -|-  1.  One  step  of  Haar 
decomposition  transforms  the  input  into  averages  Sj^k  and  details  (differences)  dj^k  at 
the  scale  j.  The  corresponding  equations  are: 


'Sj+l,2fc-|-l  ~  Sj+i^2k 

(6.10) 

Sj+l,2k  +  ’Sj+l,2fc+l 

2 

7 

Sj+l,2k  + 

(6.11) 

Thus  the  dual  lifting  step  of  the  second  generation  Haar  wavelet  decomposition  re¬ 
places  the  odd  indexed  input  values  by  the  difference  between  the  odd  and  the  even 
indexed  input  value.  Similarly,  the  primal  lifting  step  of  the  second  generation  Haar 
wavelet  decomposition  replaces  the  even  indexed  input  values  by  the  mean  of  the  odd 
and  the  even  indexed  input  value.  The  rectangular  boxes  linking  the  odd  and  the 
even  branches  in  Figure  C-50  stand  for  the  convolution  operations.  In  the  dual  lifting 
step  (Predict),  the  even  indexed  coefficients  are  convolved  with  some  sequence  (P) 
and  the  result  is  subtracted  from  the  odd  indexed  coefficients.  The  primal  lifting  step 
(Update)  convolves  the  resulting  difference  with  another  sequence  (U)  and  add  the 
result  to  the  even  indexed  coefficients.  The  convolution  operators  P  and  U  for  the 
Haar  transform  are  1  and  1/2  respectively. 


The  second  generation  wavelet  transformation  based  on  the  Haar  decomposition  is 
utilized  to  obtain  scales  in  the  chirp  signal.  A  Data  set  was  generated  at  N  [N  =  1025) 


time  instances  for  the  time  domain  (0  <  t  <  T  =  256)  using  the  function: 

3sin{u(t)t) 


where, 
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The  instantaneous  time  scale  corresponding  to  the  above  data  is  given  by: 
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32 


(6.13) 


Figure  C-51  shows  the  contour  plot  of  the  scale  coefficients  varying  with  time.  The 
scale  coefficients  are  obtained  by  the  SGW  transformation  based  on  the  Haar  decom¬ 
position  of  the  chirp  signal  given  by  Equation  6.12.  The  true  instantaneous  scales 
(Equation  6.13)  are  plotted  in  the  same  hgure  to  study  the  behavior  of  the  scale 
coefficients  obtained  using  the  SGW.  There  is  a  good  comparison  between  the  true 
instantaneous  scales  and  the  estimated  scales,  but  a  dispersion  of  the  scale  coeffi¬ 
cients  is  also  observed.  Haar  wavelets  are  an  example  of  the  simplest  wavelet,  with 
low  resolution.  However,  the  above  illustration  using  chirp  signals  shows  that  SGW 
can  potentially  be  a  powerful  tool  for  scale  estimation  in  the  ocean. 

Adaptive  scale  estimation  in  the  ocean  using  second  generation  wavelets  based 
on  other  admissible  wavelets  like  the  Mexican  Hat  or  Morkel  wavelets  is  a  subject  of 
further  research  study.  The  primal  and  dual  lifting  steps  needs  to  be  appropriately 
tuned  for  using  these  wavelets.  A  goal  will  be  to  reduce  the  dispersion  of  scale 
coefficients  for  accurate  adaptive  scale  estimation  for  oceans. 

Among  the  three  methods  discussed  above,  the  method  based  on  the  structure 
function  and  the  method  based  on  the  second  generation  wavelets  look  more  promis¬ 
ing,  since  these  methods  are  valid  for  the  non-dyadic  and  irregularly  spaced  ocean 
dataset. 
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Chapter  7 


Summary  and  Conclusions 


Our  research  consisted  of  three  related  investigations:  a.  new  methodologies  for 
the  mapping  and  dynamical  inference  of  ocean  fields  from  irregular  data  in  complex 
multiply-connected  domains,  b.  computational  studies  of  properties  of  the  new  map¬ 
ping  schemes,  and  c.  the  adaptive  estimation  of  spatial  and  temporal  scales.  Results, 
hndings  and  future  work  are  summarized  next. 

New  methods  for  efficient  field  mapping,  i.e.  Objective  Analysis,  in  complex 
coastal  regions  were  researched,  implemented  and  utilized.  These  new  OA  methods, 
which  satisfy  coastline  constraints  (e.g.  there  is  no  direct  relationship  across  land- 
forms),  are  based  on  estimating  the  length  of  the  optimal  path  using  either  the  Level 
Set  Method  (LSM)  or  the  Fast  Marching  Method  (FMM).  These  novel  methods  were 
applied  and  stndied  in  complex  domains  of  the  Philippines  Archipelago  and  Dabob 
Bay  using  realistic  datasets  to  obtain  held  estimates  snch  as  temperatnre,  salinity  and 
biology  (chlorophyll).  Results  were  compared  with  those  of  a  standard  OA  scheme 
(nsing  across-landforms  Enclidean  distance  in  the  analytical  correlation  fnnction)  and 
of  OA  schemes  based  on  the  nse  of  stochastically  forced  differential  eqnations  (SDE), 
inclnding  that  proposed  by  Lynch  and  McGillicnddy  (2001).  We  have  shown  that  onr 
new  FMM-based  scheme  is  compntationally  inexpensive  compared  to  onr  LSM-based 
scheme  and  the  SDE  approach.  Onr  illustrations  and  stndies  show  that  held  maps 
obtained  using  onr  FMM-based  scheme  do  not  reqnire  postprocessing  (smoothing)  of 
helds  e.g.  they  are  devoid  of  any  spnrions  hydrographic  held  gradients  which  are  nn- 
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acceptable  for  flow  computation.  We  have  also  shown  that  the  use  of  FMM  is  the  most 
appropriate  method  for  the  optimal  distance  estimation  among  the  distance  estima¬ 
tion  methodologies  like  Dijkstra’s  optimization  algorithm  and  the  classic  Bresenham 
line  algorithm.  The  optimal  distance  computed  using  Dijkstra’s  algorithm  is  com¬ 
putationally  expensive  and  inaccurate.  Apart  from  being  computationally  expensive, 
the  optimal  distance  computed  using  the  Bresenham  line  algorithm  is  discontinuous. 
This  results  in  the  formation  of  fronts  with  high  field  gradients.  Such  high  gradient 
fronts  do  not  occur  when  our  FMM-based  scheme  is  utilized.  We  have  also  utilized 
our  new  FMM-based  OA  scheme  for  incorporating  non-homogeneous  dynamical  ef¬ 
fects  by  appropriately  modifying  the  scalar  speed  function  in  the  Eikonal  equation. 
In  particular,  we  have  used  a  bathymetry-dependent  scalar  speed  function  to  include 
bathymetric  effects  at  lower  depth  levels.  We  also  proposed  the  use  of  the  smallest 
length  scale  on  the  optimal  path  to  include  the  non-homogeneous  effects  due  to  the 
existence  of  fronts  in  a  continental  shelf.  In  the  future,  analogous  modihcation  of  the 
scalar  speed  function  or  the  length  scale  can  be  used  to  incorporate  other  dynamical 
effects  (e.g.  conservation  of  potential  vorticity).  The  optimal  path  length  obtained 
using  our  FMM/LSM-based  scheme  can  also  be  used  to  extend  the  methodology 
proposed  by  Lermusiaux  et  ah  for  three-dimensional,  multivariate  and  multi-scale 
spatial  mapping  of  geophysical  helds  and  their  dominant  errors  (Lermusiaux  et  ah, 
1998,  2000,  Lermusiaux,  2002)  to  complex  coastal  regions.  Three-dimensional,  multi¬ 
variate  and  multi-scale  spatial  mapping  using  our  FMM  based  scheme  is  also  a  subject 
of  further  research. 

Computational  studies  of  properties  of  the  new  mapping  schemes  were  carried 
out.  The  sequential  processing  of  observations  reduces  the  computational  cost  and 
also  helps  in  understanding  the  impact  of  individual  data.  We  have  used  the  Wiener- 
Khinchin  and  Bochner  theorem  to  obtain  the  relationship  among  parameters  of  the 
analytical  correlation  function  for  it  to  be  positive  dehnite.  Such  analysis  is  valid 
only  for  the  correlation  functions  based  on  the  Euclidean  distance  for  convex  simply- 
connected  domains.  When  the  number  of  observations  is  large,  the  Kalman  Filter 
or  the  sequential  processing  of  observations  may  have  divergence  problems  due  to 
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numerical  reasons.  Methods  to  remove  divergence  problems  arising  from  numerical 
issues  were  discussed.  It  was  found  that  the  covariance  matrix  is  no  longer  positive 
dehnite  when  the  optimal  path  length  is  computed  using  FMM.  Therefore,  the  use 
of  high  order  FMM  was  discussed  and  implemented  to  obtain  more  accurate  length 
of  shortest  sea  paths.  However,  we  found  that  the  covariance  matrix  also  becomes 
negative  due  to  the  presence  of  islands  and  other  non-convex  landforms.  Several 
approaches  to  overcome  this  issue  were  discussed.  These  include  discarding  problem¬ 
atic  data  points,  introducing  process  noise,  and  reducing  the  covariance  matrix  based 
on  the  dominant  singular  value  decomposition  (SVD).  Among  these,  we  argue  that 
introducing  process  noise  and  using  the  dominant  SVD  are  the  best  solutions. 

We  have  also  derived  and  implemented  new  FMM  based  methodology  for  the  esti¬ 
mation  of  absolute  velocity  under  geostrophic  balance  in  complex  multiply-connected 
domains.  FMM  is  used  for  the  computation  of  the  minimum  vertical  area  between  all 
pairs  of  islands.  The  minimum  area  is  required  for  obtaining  the  transport  streamfunc- 
tion  which  minimizes  the  inter-island  transport  and  produces  a  smooth  velocity  flow 
held.  The  transport  streamfunction  can  then  be  utilized  to  estimate  the  geostrophic 
how  velocity  from  the  temperature  and  salinity  held  maps  alone.  We  have  illustrated 
this  method  by  applying  it  in  a  subdomain  of  the  Philippines  Archipelago. 

The  estimation  of  spatial-temporal  scales  from  irregular  ocean  held  measurements 
would  be  potentially  a  signihcant  advance  to  the  ocean  community  for  better  under¬ 
standing  and  sampling  of  ocean  processes.  This  is  a  challenging  problem  due  to  the 
multi-scale,  turbulent  and/or  intermittent  ocean  dynamics  and  due  to  the  irregular 
properties  of  the  data.  Three  new  methods  for  adaptive  spatial-temporal  scale  estima¬ 
tion  were  proposed  and  implemented.  These  methods  are  based  on  using  the  structure 
function,  short  term  Fourier  transforms  and  second  generation  wavelet  analysis.  To 
our  knowledge,  this  is  the  hrst  time  that  adaptive  methods  for  the  spatial-temporal 
scale  estimation  are  proposed  and  used  in  ocean  studies.  The  application  of  our  new 
scale  estimation  schemes  based  on  the  structure  function  method  was  illustrated  using 
Melville  (Summer  2007)  exploratory  cruise  data.  We  also  illustrated  the  application 
of  short  term  Fourier  transforms  and  second  generation  wavelet  analysis  using  chirp 


93 


signal  data.  Second  generation  wavelet  analysis  for  adaptive  spatial-time  scale  esti¬ 
mation  from  the  irregularly  spaced  ocean  data  is  shown  to  be  a  promising  technique 
and  will  be  a  subject  of  further  research. 
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Appendix  A 


Proof  of  Algorithm  for  Sequential 
Processing  of  Observations 


The  equivalence  of  Equations  5.1  —  5.3  and  Equations  5.7  —  5.9,  as  proven  by 
Parrish  and  Cohn  (1985),  is  given  here.  The  proof  will  follow  from  an  inductive 
argument  after  establishing  the  equivalence  for  J=2.  Here  we  denote  the  background 
held  and  covariance  by  and  Cor(x,x),  respectively  and  the  estimated  held  and  its 
error  covariance  by  and  respectively.  Eqnation  5.1  can  be  rearranged  to 

the  form; 


K  =  {I-  KH)Cor{x,  x)H^R-\  (A.l) 

Using  Eqnation  5.2,  Equation  A.l  can  be  written  as: 

K  =  (A.2) 

Using  Eqnation  A.2,  Equation  5.3  becomes 

^jJO^  =  4>  +  pO^H^R-\d-Hi>).  (A.3) 
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According  to  Equations  5.4  and  5.5, 


Using  A. 4,  Equation  A. 3  can  be  written  as 

^  H^R^\d2  -  H2^/J)]  (A.5) 

The  Woodbury  formula  establishes  the  equivalence  of  equations  in  A. 6  for  arbitrary 
conformable  matrices  A,  B,  C,  R,  provided  that  the  indicated  inverse  exists. 

C  =  B-  BA^{ABA^  +  R)-^AB  ^  =  5“^  +  A^R-^A  (A.6) 

The  Equations  5.1  and  5.2  lead  to 

pOA  ^  Cor{x,  x)  —  Cor{x,  x)H^{HCor{x,  x)H^  +  R)~^HCor{x,  x)  (A. 7) 

Using  A.6,  the  equivalent  statement  of  A. 7  will  be: 

(pOA^-i  _  ^Cor{x,  x))-^  +  H^R-^H  (A.8) 

The  a-posterior  error  covariance  is  obtained  by  substituting  A. 4  in  Equation  A.8 

^pOA)-i  ^  {^Cor{x,x))-^  +  +  H^R2^H2  (A.9) 

To  show  the  equivalence  of  Equations  5.1  —  5.3  and  Equations  5.7  —  5.9,  it  needs 
to  be  shown  that  A.5  and  A.9  holds  for  the  sequential  processing  of  observations 
(Equations  5.7  —  5.9).  Apply  the  Woodbury  formula  to  equations  5.7  and  5.8  and 
get 

{Cor{x,x)i)~^  =  {Cor{x,x))~^  +  H^R^^Hi  (A. 10) 
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(pOA)-i  ^  {Cor{x,x),)-^  +  H:^R^^H2  (A.ll) 

Equations  A. 10  and  A.ll  can  be  combined  to  eliminate  {Cor{x,x)i).  This  leads  to 
a-posterior  error  covariance  Equation  A. 9.  To  verify  that  A. 5  holds,  Equation  5.9  is 
written  as: 


t/’i  =  i/j  +  Cor{x,  x)iHi  (A. 12) 

=  'ilJi  +  P^^H^R2\d2-H2'ilJi).  (A.13) 

Eliminate  from  Equations  A.  12  and  A.13  to  obtain 

^  -  Hi^)  + 

P^^HjR^\d2-H2^).  (A.14) 

Premultiply  A.ll  by  P^^  and  postmultiply  it  by  Cor{x,x)  to  obtain 

pOA  =  (/  -  pO^H^R-^H2)Cor{x,x).  (A.  15) 


Using  Equation  A.  15  in  A.14  gives  Equation  A. 5.  This  proves  the  equivalence  of 
Equations  5.1  —  5.3  and  Equations  5.7  —  5.9  for  J=2.  The  proof  of  the  general  case 
follows  from  an  inductive  argument. 
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Appendix  B 


Nonlinear  least  squares  fitting 
method 


Given  a  non-linear  function  having  a  known  analytic  form  with  n  independent  pa¬ 
rameters  i.e  f{x;  /ii,  /i2,  h-3,  •••hn)  and  consider  the  over-determined  set  of  m  equations 
{m  >  n)\ 


yi  f  H2)  yZ)  ■  ■■yn) 

1/2  =  /(a^2;  hi) /^2, •••hn) 

Vm  hi ) /^2) /1'3)  •  •  •/I'n) • 

The  values  /ii,/i2,h3,---hn,  which  best  satisfy  the  above  system  of  equations  (by  min¬ 
imizing  the  sum  of  the  squared  residuals)  can  be  obtained  using  the  nonlinear  least 
squares  fitting  method.  The  approach  is  outlined  next: 

1.  An  initial  guess  for  /ij  is  chosen.  Dehne  the  error  {j3i)  as: 

Pi  Hi  f  {Xi  j  hi )  h2)  h3)  •  •  -hfi)  •^) 

2.  A  linearized  estimate  for  the  change  in  parameters  (dh*)  which  reduces  the  error 
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(yPi)  to  zero  is  given  by: 


(3  =  Adjj, 


A  = 


is  given  by 

df 

df 

dm 

d^ln 

{xx,^l) 

df 

df 

dm 

d^ln 

{x2,^l) 

df 

df 

dm 

d^ln 

Apply  the  matrix  to  both  sides  of  Equation  B.3  to  obtain 

A^  f3  =  [AA  A)d^, 


(B.3) 


(B.4) 


Equation  B.4  can  be  solved  for  djj,  using  the  standard  linear  algebraic  equation 
solvers  such  as  Gaussian  elimination  or  LU  decomposition.  The  offset  dfx  is 
applied  to  p. 

3.  The  process  is  iteratively  applied  till  the  offset  djj,  becomes  smaller  than  the 
desired  tolerance  level. 

This  is  the  nonlinear  least  squares  htting  method,  which  is  utilized  in  the  adaptive 
scale  estimation  using  the  structure  function  method. 
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Appendix  C 


Figures 
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Figure  C-1:  MSEAS  2-stage  Objective  Analysis 
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Figure  C-2:  Comparison  between  FMM  (Left)  and  the  Dijkstras  (Network  Flow) 
algorithm  (order  =  2)  (Right)  for  optimal  path  planning  (Takei,  2006) 
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Figure  C-3:  Illustration  of  the  Bresenham  Line  Algorithm 


in  situ  temperature  (°  C  )  at  0.0m 
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Figure  C-4:  Temperature  data  in  Dabob  Bay 


104 


4r5i 


47°  48 


47°  45 


47°  42 


47”  39 


122°54  122'51  122°48  122°45 


13-11 

13.02 

12.93 

12.84 

12.75 
12.66 

12.57 

12.48 
12.39 
12-30 
12.21 
12.12 
12.03 
11  94 

11.85 

11.76 
11.67 

11.58 

11.49 


122’54  122”51  122=46  122=45 


47=51 


47°  45 


47°  42 


4r39 


29-27 

29.25 

29.22 

29.20 

29.17 

29.15 

29.12 

29.10 

29.07 

29.05 

29.02 

29.00 

28.97 


28.90 

28.87 

28.85 

28.62 


122°54  122°51  122=48  122°45 


47=51 


47=48 


47°  45 


47°  42 


47°  39 


13.11 
13.02 
12.93 

12.84 
12-75 
12.66 

12.57 

12.48 
12.39 
12.30 
12.21 

12.12 
12.03 
11  94 

11.85 
11  76 
11.67 

11.58 

11.49 


122=54  122°51  122°48  122°45 


47°51 


47°  48 


4r45 


4r42 


4r39 


29.27 

29.25 

29.22 

29.20 

29.17 

29.15 

29.12 

29.10 

29.07 

29.05 

29.02 

29.00 

28.97 

28.95 

28.92 

28.90 

28.87 

28.85 

28.82 


Figure  C-5:  Temperature  (°C)  (Left)  and  Salinity  (PSU)  (Right)  OA  Fields  in  Dabob 
Bay  from  the  optimal  path  length  computed  using:  (Top)  Bresenham  Algorithm; 
(Middle)  Averaged  Bresenham  Approach;  (Bottom)  Fast  Marching  Method 
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World  Ocean  Atlas  2005  Climatology 
in  situ  temperature  (°C  )  at  0.0m 
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Figure  C-6:  World  Ocean  Atlas  2005  Climatology  in  situ  temperature  (°C)  at  0.0m 
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Figure  C-7:  Temperature  (°C)  OA  Fields  obtained  using  the  Level  Set  Method  at 
Level:  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left)  200m;  (Middle  -  Right) 
450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-8:  Temperature  (°C)  OA  Fields  obtained  using  the  Fast  Marching  Method 
at  Level:  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left)  200m;  (Middle  -  Right) 
450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-9:  Temperature  (°C)  OA  Fields  (Standard  OA  without  taking  islands  into 
account)  at  Level;  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left)  200m;  (Middle 
-  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-10:  Temperature  (°C)  OA  Fields  using  the  SDE  approach  (representing 
covariance  by  Helmholtz  equation)  at  Level:  (Top  -  Left)  Om;  (Top  -  Right)  40m; 
(Middle  -  Left)  200m;  (Middle  -  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  - 
Right)  3000m 
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Figure  C-11:  Temperature  (°C)  OA  Fields  using  the  SDE  approach  (representing  held 
by  Helmholtz  equation)  at  Level;  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left) 
200m;  (Middle  -  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-12:  Salinity  (PSU)  OA  Fields  obtained  using  the  Level  Set  Method  at  Level: 
(Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left)  200m;  (Middle  -  Right)  450m; 
(Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-13:  Salinity  (PSU)  OA  Fields  obtained  using  the  Fast  Marching  Method  at 
Level:  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left)  200m;  (Middle  -  Right) 
450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-14:  Salinity  (PSU)  OA  Fields  (Standard  OA  without  taking  islands  into 
account)  at  Level;  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left)  200m;  (Middle 
-  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-15:  Salinity  (PSU)  OA  Fields  using  the  SDE  approach  (representing  covari¬ 
ance  by  Helmholtz  equation)  at  Level:  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  - 
Left)  200m;  (Middle  -  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-16:  Salinity  (PSU)  OA  Fields  using  the  SDE  approach  (representing  held 
by  Helmholtz  equation)  at  Level;  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left) 
200m;  (Middle  -  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-17:  Comparison  of  the  Temperature  (°C)  field  at  Level  =  1000m  by  using: 
(Top  -  Left)  Level  Set  Method;  (Top  -  Right)  Fast  Marching  Method;  (Bottom  -  Left) 
SDE  Approach;  (Bottom  -  Right)  Standard  OA 
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Figure  C-18:  Comparison  of  the  Salinity  (PSU)  field  at  Level  =  1000m  by  using: 
(Top  -  Left)  Level  Set  Method;  (Top  -  Right)  Fast  Marching  Method;  (Bottom  -  Left) 
SDE  Approach;  (Bottom  -  Right)  Standard  OA 
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Figure  C-19:  Difference  between  Temperature  (°C)  field  at  Level  =  1000m  obtained 
using  Fast  Marching  Method  and  using:  (Top  -  Left)  Level  Set  Method;  (Top  -  Right) 
Standard  OA;  (Bottom  -  Left)  SDE  Approach  (representing  covariance  by  Helmholtz 
equation);  (Bottom  -  Right)  SDE  (representing  held  by  Helmholtz  equation) 
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Figure  C-20:  Difference  between  Salinity  (PSU)  field  at  Level  =  1000m  obtained 
using  Fast  Marching  Method  and  using:  (Top  -  Left)  Level  Set  Method;  (Top  -  Right) 
Standard  OA;  (Bottom  -  Left)  SDE  Approach  (representing  covariance  by  Helmholtz 
equation);  (Bottom  -  Right)  SDE  (representing  held  by  Helmholtz  equation) 
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Figure  C-21:  Melville  exploratory  Cruise  +  GTSPP  +  HB2  Climatology  (Summer 
2007)  in  situ  temperature  (°C)  at  0.0m 
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Figure  C-22:  Temperature  ("C)  and  Salinity  (PSU)  Field  Maps  (Melville  exploratory 
Cruise  +  GTSPP  +  HB2  Climatology  (Summer  2007)  data) 
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Figure  C-23:  Melville  exploratory  cruise  and  glider  data  (Summer  2007)  in  Philippines 
Archipelago 
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Figure  C-24:  Temperature  (°C)  OA  Fields  using  the  FMM  with  Melville  exploratory 
cruise  and  glider  data  (Summer  2007)  at  Level:  (Top  -  Left)  Om;  (Top  -  Right)  40m; 
(Middle  -  Left)  200m;  (Middle  -  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  - 
Right)  3000m 
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Figure  C-25:  Salinity  (PSU)  OA  Fields  using  the  FMM  with  Melville  exploratory 
cruise  and  glider  data  (Summer  2007)  at  Level:  (Top  -  Left)  Om;  (Top  -  Right)  40m; 
(Middle  -  Left)  200m;  (Middle  -  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  - 
Right)  3000m 
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PHILEX  -  Joint  Cruise  -  Melville  Data 
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Figure  C-26:  Philippines  Archipelago  -  Melville  joint  cruise  Data  (Winter  2008) 
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Figure  C-27:  Temperature  (°C)  OA  Fields  using  the  FMM  with  Melville  joint  cruise 
data  (Winter  2008)  at  Level;  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left) 
200m;  (Middle  -  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-28:  Salinity  (PSU)  OA  Fields  using  the  FMM  with  Melville  joint  cruise  data 
(Winter  2008)  at  Level;  (Top  -  Left)  Om;  (Top  -  Right)  40m;  (Middle  -  Left)  200m; 
(Middle  -  Right)  450m;  (Bottom  -  Left)  1000m;  (Bottom  -  Right)  3000m 
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Figure  C-29:  Biology  (chlorophyll)  data  in  Philippines  Archipelago 
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Figure  C-30:  Chlorophyll  (/iinol/Kg)  OA  Fields  using  the  FMM  at  Level:  (Top  - 
Left)  Om;  (Top  -  Right)  10m;  (Middle  -  Left)  50m;  (Middle  -  Right)  100m;  (Bottom 
-  Left)  160m;  (Bottom  -  Right)  200m 
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Figure  C-31:  Velocity  estimation  under  geostrophic  balance  (weight  functions  based 
on  the  minimum  vertical  area)  from  held  maps  (WOA05)  obtained  using  the  FMM: 
(Top  -  Left)  Streamfunction,  Velocity  at  depths:  (Top  -  Right)  Om;  (Bottom  -  Left) 
100m;  (Bottom  -  Right)  1000m 
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Figure  C-32:  Velocity  estimation  under  geostrophic  balance  (weight  functions  based 
on  the  minimum  vertical  area)  from  held  maps  (WOA05)  obtained  using  the  SDE 
(Helmholtz  equation)  for  held:  (Top  -  Left)  Streamfunction,  Velocity  at  depths:  (Top 
-  Right)  Om;  (Bottom  -  Left)  100m;  (Bottom  -  Right)  1000m 
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Figure  C-33:  Velocity  estimation  under  geostrophic  balance  (weight  functions  based 
on  the  minimum  distance)  from  held  maps  (WOA05)  obtained  using  the  FMM:  (Top 
-  Left)  Streamfunction,  Velocity  at  depths:  (Top  -  Right)  Om;  (Bottom  -  Left)  100m; 
(Bottom  -  Right)  1000m 
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Figure  C-34:  Velocity  estimation  under  geostrophic  balance  (weight  functions  based 
on  the  minimum  distance)  from  held  maps  (WOA05)  obtained  using  the  SDE 
(Helmholtz  equation)  for  held:  (Top  -  Left)  Streamfunction,  Velocity  at  depths:  (Top 
-  Right)  Om;  (Bottom  -  Left)  100m;  (Bottom  -  Right)  1000m 
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Figure  C-35:  Cartoon  illustrating  non-homogeneous  scales  in  a  continental  shelf  with 
a  front  having  small  scales  in  y-direction  separating  regions  having  large  scales  in 
y-direction 
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Figure  C-36:  Inclusion  of  bathymetry  effects  by  modifying  scalar  speed  function  in 
Philippines  Archipelago  at  a  depth  of  450m:  (Top  -  Left)  Regular  scalar  speed  function 
(FI)  with  0  for  landforms  and  1  for  ocean;  (Top  -  Right)  Modihed  scalar  speed  function 
(F2)  given  by  bathymetry  /  OA  level;  (Middle  -  Left)  Temperature  (°C)  held  map 
based  on  FI;  (Middle  -  Right)  Temperature  (°C)  held  map  based  on  F2;  (Bottom  - 
Left)  Salinity  (PSU)  held  map  based  on  FI;  (Bottom  -  Right)  Salinity  (PSU)  held 
map  based  on  F2 
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Figure  C-37:  World  Ocean  Atlas  2005  (Spliced  February  and  Winter  Climatology)  in 
situ  temperature  ('’C)  at  0.0m 
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Figure  C-38:  Temperature  ("C)  OA  Fields  using  the  Fast  Marching  Method  at  the 
surface  (Om)  using  the  following  scales:  (Left)  Lq  =  540iFm,  =  ISOKm;  (Right) 

Lq  =  lOSOKm,  Le  =  360iFm 
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Figure  C-39:  Comparison  between  distance  and  correlation  {exp{^);L  =  10)  ob¬ 
tained  using  the  first  order  FMM  and  the  Euclidean  distance  on  a  30-by-30  domain 
without  islands.  (Top  -  Left)  Difference  in  distance;  (Top  -  Right)  Normalized  differ¬ 
ence  in  distance;  (Bottom  -  Left)  Difference  in  correlation;  (Bottom  -  Right)  Normal¬ 
ized  difference  in  correlation 
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Figure  C-40:  Comparison  between  distance  and  correlation  {exp{^);L  =  10)  ob¬ 
tained  using  the  higher  (second)  order  FMM  and  the  Euclidean  distance  on  a  30-by-30 
domain  without  islands.  (Top  -  Left)  Difference  in  distance;  (Top  -  Right)  Normal¬ 
ized  difference  in  distance;  (Bottom  -  Left)  Difference  in  correlation;  (Bottom  -  Right) 
Normalized  difference  in  correlation 
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Figure  C-41:  Temperature  (°C)  OA  Fields  at  the  surface  (Om)  (scales  Lq  =  lOSOKm, 
Lf,  =  360iFm)  using:  (Left)  First  Order  FMM;  (Right)  Higher  order  FMM 
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Figure  C-42:  Example  of  an  idealized  (multiply-connected)  domain  having  an  island 
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Figure  C-43:  Example  illustrating  the  mapping  and  stretching  in  a  domain  with  a 
circular  island.  Multiply-connected  domain  in  the  polar  (r,^)  coordinates  is  mapped 
to  a  simply-connected  rectangular  domain  in  the  curvilinear  (C,h)  coordinates. 
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Figure  C-44:  Temperature  (°C)  OA  Fields  at  the  surface  (Om)  (scales  Lq  =  lOSOKm, 
Le  =  360Km)  using:  (Top  -  Left)  FMM;  (Top  -  Right)  FMM  and  by  discarding 
the  problematic  data;  (Bottom  -  Left)  FMM  and  by  introducing  process  noise;  (Bot¬ 
tom  -  Right)  FMM  and  by  dominant  singular  value  decomposition  (SVD)  of  a-priori 
covariance 
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Figure  C-45:  Scale  estimates  for  Philippines  Archipelago  from  Melville  exploratory 
cruise  data  (Summer  2007)  on  26*^  June  2007  (Left)  and  15*^  July  2007  (Right)  using 
structure  function 
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Figure  C-46:  Small  and  large  length  scales  (Top)  and  time  scales  (Bottom)  from  the 
adaptive  learning  algorithm  (learning  rate  =  0.1)  based  on  STFT  applied  on  the  chirp 
signal  data 
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Figure  C-47:  Small  and  large  length  scales  (Top)  and  time  scales  (Bottom)  from  the 
adaptive  learning  algorithm  (learning  rate  =  0.2)  based  on  STFT  applied  on  the  chirp 
signal  data 
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Figure  C-48:  Small  and  large  length  scales  (Top)  and  time  scales  (Bottom)  from  the 
adaptive  learning  algorithm  (learning  rate  =  0.1)  based  on  STFT  applied  on  Data 
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Figure  C-49:  Haar  Wavelet  (Left)  and  Mexican  Hat  Wavelet  (Right)  (Jansen  and 
Oonincx,  2005) 
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Figure  C-50:  Block  diagram  of  second  generation  wavelet  transform  (Jansen  and 
Oonincx,  2005) 
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Figure  C-51:  Scale  coefficients  for  idealized  chirp  signal  using  second  generation 
wavelets  (Lifting  scheme  corresponds  to  Haar  Wavelet).  Black  line  is  the  plot  of 
instantaneous  scale  of  the  chirp  signal  varying  with  time. 
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